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Frequency hopped signals are widely used in various communication applica- 
tions for their inherent security features. The demand, by civilian and governmental 
agencies, to intercept communication signals is increasing. The interception task can 
be summarized by detecting the signal’s presence in additive noise, classifying the 
modulation type, estimating the control parameters, decoding the data, and decrypt- 
ing the information content. This work addresses the merging of wavelet and corre- 
lation concepts to detect, classify and estimate the parameters of frequency hopped 
signals. 

We address the interception problem in two ways. The first approach is based 
on a visual inspection of the wavelet surfaces generated from the instantaneous cor- 
relation function of the communication signal and leads to hop start/stop times es- 
timates. In the second approach, we apply an energy-based processing scheme to 
estimate the hop start and stop times, the hop-scale pattern, and the hop frequency. 

Results show that frequency hopped signals can be identified at an SNR of 3 
dB or 6 dB using visual inspection or an automated scheme, respectively. 
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I. 



INTRODUCTION 



A. PROBLEM STATEMENT AND OVERVIEW 

The demand to be able to intercept communication signals is increasing, as 
interception is a basic investigation tool for civilian, intelligence and military author- 
ities. The task of intercepting a communication signal can be summarized by detect- 
ing the signal’s presence in additive noise, classifying the modulation type, estimating 
the control parameters required for reception, decoding the data, and decrypting the 
information content. Sometimes, the procedure is stopped at the detection (identifi- 
cation) step such as in geo-location or spectrum monitoring. Intelligence applications 
aim at the final step which is to extract the information content. 

Communication signals can utilize a great diversity of digital modulation tech- 
niques. Among the digital techniques, the spread spectrum signal (SS) is widely used. 
Frequency hopping signals (FH) are a subset of SS, and are used primarily in this 
dissertation. In the open literature, many signal processing tools are available for the 
interception task, and among those are the correlation analysis and wavelet analysis. 
Each of these tools has been used independently in the signal processing and commu- 
nication area. In this work we will address the merging of the wavelet and correlation 
concept to detect, classify and estimate signal parameters. 

The FH signal is a non-stationary process. Hence a specific type of two- 
dimensional correlation function, called the instantaneous correlation function, is se- 
lected, as it accommodates the time- varying nature of the signal of interest. Ap- 
plication of wavelet analysis to correlation functions is a new area and is still in 
the exploratory stage. The automation of the interception and exploitation of dig- 
ital communication signals is the final goal. This work addresses the interception, 
identification, classification, and parameter extraction of FH signals. 

We start by investigating the wavelet transform of different types of the cor- 
relation functions. Then, we select the instantaneous correlation function (ICF) as 
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the candidate correlation representation. We derive the instantaneous correlation 
function of the FH signal and analyze its wavelet transform obtained by using the 
Morlet wavelet. The wavelet transform of the instantaneous correlation function is 
a 3-dimensional (3-D) surface. We partition the 3-D surface into a number of 2-D 
surfaces, where each corresponds to one wavelet scale (called the wavelet surface). 

We address the interception problem in two approaches. First, we visually 
inspect the 2-D wavelet surfaces to identify and classify the structure of the FH 
signal and obtain an estimate for the hop time interval. Second, we apply a proposed 
processing scheme to estimate the hop start/stop time, the hop-scale pattern, and the 
hop frequency. The extraction of the hop start/stop time is addressed using an edge 
detection approach by applying a compass operator, a technique which is well known 
in the image processing area. The hop-scale pattern is obtained by applying an energy 
analysis. The energy analysis assigns a scale index (called the proper scale) to each 
hop. The proper scale of each hop is that scale which has the greatest energy content. 
The sequence of proper scales, representing the hop sequence, is called the hop-scale 
pattern. The frequency of each hop can be extracted from the wavelet surface at 
the proper scale. Thus, the identification and classification of the FH signal may be 
accomplished as follows: 

1. Based on the hop-scale pattern: 

If two or more wavelet scales are applied, the hop-scale pattern of the FH 
signal will be different from the hop-scale patterns of other digital modulation 
signals. 

2. Based on the frequency diversity: 

If all hop frequencies reside in one scale then the FH signal will have different 
frequency components as a function of the hop intervals. 

This dissertation is organized into eight chapters. In Chapter I, we introduce 
the problem and review related work with a focus on the interception of FH signals. 
In Chapter II, we introduce the most significant signal analysis tools. In Chapter III, 
we derive the wavelet response to correlation functions. In Chapter IV, we derive the 
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ICF for the FH signal. In Chapter V, we perform the wavelet analysis of the ICF of 
the FH signal. In Chapter VI, we analyze the processing scheme. In Chapter VII, 
we provide the simulation results. Finally, Chapter VIII contains conclusions and 
recommendations for future work. 

B. INTERCEPTION OF DIGITAL COMMUNICATION 
SIGNALS 
1. Introduction 

Interception of communication signals is of interest to a wide range of appli- 
cations in surveillance, intelligence, reconnaissance, geo-location, spectral monitoring 
and jamming [Ref. 1]. 

Digital communication systems can use a large diversity of modulation tech- 
niques. Examples of these techniques are: ASK, BPSK, BFSK, QAM, MPSK and 
MFSK. Interception of digital communication signals consists of detection (identifi- 
cation), classification (feature extraction), parameter estimation, decoding, and de- 
crypting. Through out this dissertation we will use the following definitions: 

1. Detection or identification is the process of discriminating between noise and 
a digitally modulated signal. 

2. Classification (feature extraction) is the process of discriminating between dif- 
ferent modulation types, by obtaining significant parameters or signal constel- 
lation which uniquely specify the modulation type. 

3. Parameter estimation is the process of extracting the control parameters of 
the signal. 

4. Decoding is the process of detecting and correcting the errors of the commu- 
nication channel. 

5. Decrypting is the process of restoring the original information content of the 
transmitted messages. 
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2. Interception Tools 

A large number of open-literature publications address the interception of 
digital communication signals. The different signal processing tools used for signal 
interception may be categorized as follows: 

1. Basic tools: Spectral analysis, correlation analysis, and parametric modeling. 

2. Linear tools: Linear transforms including the wavelet transform. 

3. Non-linear tools: Higher-order spectra, spectral correlation, and cyclic-feature 
processing. 

4. Other tools: Eigen-analysis, singular-value decomposition, and stochastic res- 
onance. 

There are two different approaches to address the interception. The first approach 
addresses the interception as a decision-making problem applied to a set of alternative 
hypotheses. This requires statistical decision theory and hence a statistical description 
of the alternatives. The second approach handles the interception as a statistical 
pattern-recognition problem. 

3. Digital Signal Modulation 

Different tools can be used for intercepting digital modulations. For ASK, 
QAM, BPSK, QPSK, and FSK, applications are addressed in references [Ref. 2] to 
[Ref. 14], [Ref. 17] and [Ref. 32]. Since the FSK modulation is related to the FH 
signal, the interception tools for FSK are briefly introduced as follows: 

Spectral correlation and visual inspection is applied for identification of FSK signals 
in [Ref. 2]. Higher-order moments are applied for FSK signals in [Ref. 3]. Wavelet 
analysis of the FSK signal is applied for timing extraction in [Ref. 5, 17] where it is 
shown that the frequency transitions of the FSK signal are related to inflection points 
contained in the scalogram. 
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4. Spread Spectrum Signals 

Interception of spread spectrum signals is addressed in [Ref. 18] to [Ref. 35]. 
There are five conventional techniques to intercept frequency hopped (FH) signals 
[Ref. 18, 19]. These are: 

1. Wideband energy detectors (vsdde band radiometer), 

2. Optimum-multichannel FH pulse-matched filters, 

3. Maximum-channel filter-bank combiners, 

4. Optimum partial-band FH pulse-matched detectors, and 

5. Partial-band filter-bank combiners. 

These techniques differ mainly in two points. The first difference is in the bandwidth 
of the interception filter(s) relative to the bandwidth of the FH signal. The second 
difference is in the number of parallel channels relative to the number of hopping 
frequencies. Consequently, they differ in the required minimum signal-to-noise ratio 
(SNR) for acceptable performance. 

Time or frequency uncertainty can be minimized by using overlapping tech- 
niques [Ref. 20]. A likelihood-ratio-test detector and a channelized receiver are widely 
accepted as the optimum system for detecting FH signals [Ref. 20] to [Ref. 24]. The 
maximum-likelihood detector (ML), using a bank of correlators, is addressed in [Ref. 
22]. It assumes prior knowledge of the hopping times and frequencies, and one corre- 
lator is designated to each primary frequency region. Results show that the coherent 
ML scheme gives a probability of detection of 0.5 at an SNR of 4.5 dB for a probabil- 
ity of false alarm of 10“®. For the noncoherent ML scheme, this performance requires 
an SNR of about 5.9 dB. 

A generalized likelihood ratio test for a multiple hop observation interval is 
addressed in [Ref. 23]. It assumes prior knowledge of hop frequencies, hop rate, 
and hop times. Using an observation interval of 10^ hops and a probability of false 
alarm of 10“^, a probability of miss of 10“^ is achieved at an SNR of 8 dB. The same 
performance requires an SNR of 5 dB using an observation interval of 10® hops. 



5 



Applying wavelet analysis to the interception of FH signals is addressed in 
[Ref. 35]. FH classification is based on locating transition spikes due to frequency 
transitions. The authors suggest using the STFT for the hop frequency estimation 
instead of the wavelet transform because wavelets have varying frequency resolution 
and the FH signal has a wide bandwidth. 
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II. 



FOURIER ANALYSIS AND WAVELETS 



Signal analysis (expansion, decomposition, or transformation) is a method 
used to represent time signals as a linear combination of elementary building blocks 
(or elementary basis functions). This representation is vital to the area of signal 
processing. Each expansion is defined by its basis functions, or equivalently, its 
basic building blocks. There are numerous linear expansion methods. Well-known 
examples are the Shannon expansion, the Karhunen-Loeve expansion, the Gram- 
Schmidt expansion, the eigen-analysis, and the most popular, the Fourier analysis. 
In this chapter, we summarize the evolution of the analysis tools in three steps. First, 
we introduce the Fourier analysis as used for stationary signals. Then, we introduce 
time-frequency distributions, which can be used to represent time-varying signals. 
Finally, we introduce the concept of the wavelet analysis to represent time-varying 
signals over a time-scale plane. 

A. FOURIER ANALYSIS 

The major reference for this section is [Ref. 36]. The Fourier transform (FT), 
also called Fourier analysis, is the most popular signal decomposition (expansion). 
The Fourier transform is used to represent a stationary process by decomposing it into 
sinusoidal or complex exponential components. Here, a stationary process is defined 
as a process whose spectrum is not varying with time. Hence, Fourier techniques work 
well and allow successful frequency localization of the spectral components. When 
a non-stationary process is present, its time- varying spectrum requires an additional 
dimension (i.e., time). The major problem with the classical Fourier analysis is that 
it uses an infinitely long sinusoidal or complex exponential basis of functions. This 
makes time localization impossible. 

By classical Fourier analysis we mean the Fourier series (FS), the continuous- 
time Fourier transform (FT), the discrete-time Fourier series (DFS), and the discrete- 
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time Fourier transform (DTFT). 

The discrete Fourier transform (DFT), or its fast implementation, the fast 
Fourier transform (FFT), uses finite integration time. The Fourier transform cannot 
give a time resolution better than the integration interval. To cope with the re- 
quirement of tracking time evolution and to provide time localization, the short-time 
Fourier transform (STFT) was introduced. The STFT uses a window function which 
affects the time resolution; that is, the longer the window interval the worse is the 
time resolution. 

1. Fourier Series 

Given a periodic continuous time signal x{t), with period T, the Fourier series 
(FS) expansion is given by : 

CX) 

x{t)= ^ (II.l) 

k=-—oo 

where 

C{k) = (II.2) 

and fo = 

Note that the expansion coefficients C{k) are evaluated at integer values of k. 
Parseval’s relation for this power signal ( recall that power signals have finite average 
power and infinite energy) is given by : 

p. = i / |x(t)|Vt = f; |c(*)|^ (II.3) 

^ k=-oo 

where Px is the average power of the signal. The signal has to satisfy the Dirichlet 
conditions to have a valid Fourier series expansion. These conditions can be summa- 
rized as: the signal has a finite number of discontinuities, a finite number of maxima 
and minima, and it is absolutely integrable (over one period). 
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2. Fourier Transform 

A non-periodic continuous time signal, x(t), can be represented as 

/ oo , 

XU)e>^‘'df, (II.4) 

-00 

where 

(II.5) 

-00 

Note that both the signal and the transform are continuous functions of time and 
frequency, respectively. 

The Dirichlet conditions have to be met by the non-periodic signal with one 
modification; that is, the time support of the signal is considered instead of the period 
in the above-addressed Dirichlet conditions. 

Parseval’s relation is given by 

E. = r \x(t)fdt = r \x{f)\'‘d/ (II.6) 

J — OO 00 

where is the total energy of the signal. 

Short-Time Fourier Transform. Because of the lack of time localization, non- 
stationary signals have time-varying spectra which cannot be represented by the clas- 
sical Fourier transform. The time evolution of the time-varying spectra needs to be 
considered. Hence, there is a need for an expansion in time and frequency such as the 
short-time Fourier transform (STFT). This transform windows the signal around a 
certain time instant, performs the frequency domain analysis, and repeats the process 
at other time instants. Note that it is assumed that the windowed signal will have 
a non-time-varying spectrum (local stationarity) within the time window. Therefore, 
the STFT for a continuous signal x(t) is given by 

/ OO 

x{t)g*{t - (II.7) 

-00 

where g{t) is the window function. X{uj, 1) is the spectral description of x{t) when the 
window is centered at the time 1. When the window g[t) has a Gaussian shape, the 
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STFT is called the Gabor transform. The STFT depends on a window function with 
a fixed width (in time). Recall that the time localization and frequency localization 
of the STFT are controlled by the window width. Consequently, the STFT offers 
a fixed-time resolution and a related fixed-frequency resolution, which results in a 
uniform tiling in the time-frequency plane. 

The time and frequency resolution are governed by the uncertainty principle, 
which will be addressed shortly. For details, see [Ref. 38]. The STFT for discrete 
time signals is defined as 

OO 

X{k,n)= g(m — n)x(m)e~^^^^M ^ (H-8) 

m=oo 

where A: = 0, 1, . . . , M — 1. 

3. Discrete-Time Fourier Series 

The discrete-time Fourier series (DTFS) for a periodic discrete signal x{n) 
with period N is given by 

N-l 

x{n) = Y, (II.9) 

k=0 

and 

C{k) = , (11.10) 

where C{k) represents the expansion coefficient at the discrete frequency Wk = 2Trk^. 
Parseval’s relation is given by 

fc =0 n =0 

where P* is the average power of the infinite energy periodic signal x{n). The Dirichlet 
conditions must be met by x{n). 

4. Discrete-Time Fourier Transform 

The discrete-time Fourier transform (DTFT) for finite length non-periodic 
discrete signals is given by 

x{n) = ^ f X{u})eP‘^^dw, (11.12) 

ZTT J2Tr 
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where 



X{u)= ^ (11.13) 

k=-—oo 

Note that the discrete frequency u) ranges from — tt to tt. Parseval’s relation 
is given by 

oo 1 . 

E.= Y: kNI" = ^ / (n.14) 

n=-oo 27T J 

where Ex is the total energy of the signal. 

Discrete Fourier Transform. A finite-length non-periodic discrete signal has a 
continuous-ffequency-domain representation X{u). For digital signal processing it is 
easier to represent X{u) by its discrete-frequency samples. This leads to the discrete 
Fourier transform (DFT). The DFT is given by 

*=0 

and 

X(k) = Z x(n)Wi5”. (H-16) 

71=0 

where and n,k = 0, 1, • • • , AT — 1. The DFT requires complex 

multiplication operations and N{N—\) complex addition operations. The fast Fourier 
transform (FFT) implements the DFT with fewer multiplications. The FFT has 
computation complexity (number of multiplications) ^Ioq 2 N if N is power of 2. 

B. TIME-FREQUENCY DISTRIBUTIONS 
1. Introduction 

A time-frequency distribution (TFD) describes non-stationary signals by dis- 
playing the energy density over time and frequency simultaneously. Let the signal 
and its Fourier transform be denoted by ${t) and S{u), respectively. The following 
terms and definitions [Ref. 38] are useful in TFD discussions: 

• Energy density or instantaneous power is given by 

|s(i)p. (11.17) 
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• Total energy of the signal is given by 



J \s{t)f dt. 


(11.18) 


Mean time is given by 

{t) = J t\s{t)fdt 


(11.19) 


Time variance is given by 


(11.20) 


where a/ is an indication of the duration of the signal. 
Energy density spectrum is given by 

|5(o;)p. 


(11.21) 


It is the intensity per Hertz near the frequency u. 
Parseval’s theorem states that 

J \S{u))f du! = j |s(t)pdt. 


(11.22) 


Mean frequency is given by 

(w) = J u\S{u})\‘^ duj . 


(11.23) 


Frequency variance is given by 


(11.24) 



where is an indication of the rms bandwidth. 

• The uncertainty principle: 

It is a fundamental statement regarding the Fourier transform pairs. It governs 
the relationship between the spread of any signal in the time domain and in 
the frequency domain. This relationship is given by atcr^ > 1/2. This means 
that there is a trade-off time localization and frequency localization. 
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2. Fundamental Properties of Time-Frequency 
Distributions 

The major reference for these properties is [Ref. 38]. Let p{t,u) be the inten- 
sity of the signal s{t) at time t and frequency to. Then p{t,uj)AtAuj is the fractional 
energy in the cell AtAu at time t and frequency u. Basic properties of the TFD are 
given by: 

• Marginal property: Summing up the energy distribution for all frequencies at 
a particular time gives the instantaneous energy. Summing over all times at a 
particular frequency gives the energy density spectrum, i.e., 

J p{t,u)du> = |s(t)|^ (11.25) 

and 

J p{t,uj)dt= |S'(o;)p. (11.26) 

Note: Any joint distribution that yields the correct marginal property is con- 
sistent with the uncertainty principle. 

• Total energy property: The total energy of the distribution should be equal 
to the total energy of the signal. 

E = I Jp{t,u)cliodt = I |s(t)|^dt = I |5(w)|^da; (11.27) 

Note: If the distribution satisfies the marginal property then the total energy 
property is automatically satisfied. 

• Time-shift-invariance: The shifted signal s{t — to) will have the distribution 
p{t-to,u). 

• Frequency-shift-invariance: The modulated signal s{t)e^^*' will have p{t, cv — ^) 
as its distribution. 

• Time and frequency shift invariance: The signal s{t — to)e^^^ will have p{t — 
to, u — ^) as its distribution. 

• Linear scaling property: The signal Ssc{t) = y/as{at) will have the distribution 
Psc{tjto) = p{at,u/a), which will satisfy the marginal property of the scaled 
signal Sic(0- 

• Weak finite-support property: If s(t) has zero value outside the time interval 
{ti, t 2 ) or outside the frequency range (o;i, U 2 ), so does the distribution p{t, lo). 
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• Strong finite-support property: If s{t) has zero value at any time to or at any 
frequency coq, then p{t, u) will be zero for any time to or frequency uq- 

Note: It is impossible for the signal f{t) to be both time and frequency 
limited. Therefore, if a distribution satisfies the weak finite-support property, 
it cannot be limited to a finite region of time-frequency plane. 



3. Wigner-Ville Distribution 

The major reference of this topic is [Ref. 38, 39]. The reason the Wigner-Ville 
distribution (WVD) is discussed in detail is that we relate the WVD and the wavelet 
transform of the instantaneous correlation function. Another reason is that the WVD 
is typically used to represent a non-stationary process. 

The WVD, in terms of the signal s{t) or its spectrum S{o}), is defined by: 



= [^1 /->* (^H^) 



't -f- r' 



(11.28) 



and 






' (j) -H 6^ 



' LO — 6^ 






(11.29) 



4. Properties of the Wigner-Ville Distribution 

The WVD has the following properties: 



• The WVD is always real, even if s{t) is a complex signal. This is because 

• The WVD satisfies the following symmetry relations: 

W{t,u>) = W{t,—u) if s{t) is real, i.e., if S{u)) = S{—co), 



and 



W(t, w) = W(— t, w) if S{u) is real, i.e., if s{t) = 

• The marginal property is satisfied in both time and frequency. Hence the total 
energy property is also satisfied. 

• The time and frequency shift invariance properties are satisfied. 
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• For s{t) = the conditional spread in frequency is given by 



oi,t = - mtf = 0.5 





(11.30) 



where denotes the time derivative. Since Equation 11-30 may result in neg- 
ative values, it cannot be properly used as a measure of bandwidth. Therefore, 
while WVD can be used to compute a reliable mean frequency, it cannot be 
used to obtain reliable values for the spread of those frequencies. 

• WVD of sum of two signals: If s{t) = Si{t) + S 2 {t), then the WVD is given by 



W{t, 0 j) = Wn{t,u) + W22{t,u) + Wi2{t,oj) + W2i{t,oj), (11.31) 



where 






is the cross WVD. The cross WVD is complex-valued, but W 12 = W^i- Hence, 
W\2{t,uj) -f W2i{t,oj) is real and we can express W{t,u) as: 



W{t,uj) = -t- W22{t,uj) + 2Re [Wi 2 (t,o;)] . (11.33) 



The additional term 2Re [IVi 2 (t,w)] is often called the cross term or the inter- 
ference term. 

Note: The signal is said to be composed of more than one signal (i.e., multi- 
component signal) if it can be segmented into separate regions in the time-frequency 
plane. 

Note: The WVD at time t depends on the values of the signal at the time 
moments [t — r/2) and {t -I- r/2). Even if there is no noise at time t, the WVD will 
reflect the noise component existing at those time moments (t— r/2) and {t + t/2). 
There are other types of TFDs, all forming one class, called the general Cohen’s 
class. Cohen’s class incorporates different kinds of TFDs which differ only in the 
definition of the window function. Window functions may be chosen to optimize the 
TFD according to certain criteria. One example of this criteria might be to minimize 
cross terms. 
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C. WAVELET ANALYSIS 



Wavelet analysis is a new trend for representing non-stationary signals. It is 
widely applied for many areas of signal processing. In the following discussion we 
introduce the basic concepts of wavelet analysis. 

1. Frames 

A basis for a vector space F is a set S of vectors in V such that S is linearly 
independent and S spans V [Ref. 41]. Therefore, a basis can express any vector 
in the space F by a linear expansion. The resultant set of expansion coefficients is 
unique. Completeness and uniqueness of the expansion coefficients are directly related 
to the exact reconstruction of the original vector from the expansion coefficients. 
This property is required in signal processing applications which use the expansion 
algorithm to reconstruct the original signals. For simple computation of the expansion 
coefficients, the vectors of S are required to be orthogonal. Orthogonality will make it 
possible to compute the expansion coefficients (analysis) using simple inner products 
with the basis vectors. Reconstruction (synthesis) and decomposition (analysis) are 
done using the same set of basis vectors. 

Orthogonality of the basis vectors (or basis functions) is defined as follows: 
the set of basis vectors Vk{t) is said to be orthogonal if 

{Vk{t), Vm{t)) = ak,m^{k - m), (11.34) 

where A:, m are integers, ak,m is a constant and S is the Kronecker delta function. The 
delta function S{k — m) will be zero everywhere except where m = k it will be one. 
Since the orthogonality is not a condition for having a set of basis vectors, we may 
find a set of non-orthogonal basis vectors. For any orthogonal set of basis vectors Vk{t) 
the finite energy time signal f{t) can be expanded using the set of basis functions 
Vk{t) as: 
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where 



% = {f{t),Vk{t)), 

and is the norm of the vector Further, if we assume unit norm for 

the vectors Vk{t), the vectors become an orthonormal set, and the expansion formula 
will be 

( 11 . 36 ) 

k 

where 

(^k = {f{t),Vk{t)). 

It is worth mentioning here that due to orthogonality (or orthonormality) both 
the analysis and synthesis of the expansion are carried out using the same set of basis 
vectors. 

For a biorthogonal basis, we need two different sets of basis vectors, the original 
set Vk{t) (the basis) for the expansion and a dual set Vk{t) (dual basis) for reconstruc- 
tion (or synthesis). The sets of the basis and the dual basis are not orthogonal, but 
each vector and its dual are orthogonal. The expansion is given by 

f{t) = J^dkVk{t), (11.37) 

k 

where 

If the set S of the expansion vectors does not satisfy the linear independence 
condition stated above, then the set is called a “frame.” Expansion with frames does 
not maintain uniqueness of expansion coefficients [Ref. 37]. Also, the exactness of 
the reconstruction will be replaced with an approximate representation. For more 
details about the theory of frames see [Ref. 42, 43]. Frames do not satisfy the 
Parseval’s theorem of energy, which states that the energy in the original function 
will be distributed among the expansion coefficients without loss or gain. For frames, 
the total energy in the expansion coefficients will have two bounds A and B such 
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that: 



< EI(/W.»‘W)P < BII/IP (11-38) 

A: 

where 0 B < oo. Note that \i A = B, then the frame is called a tight frame and 

>1||/IP = E !(/(*). (II.39) 

k 

If ^ = S = 1, then the frame becomes an orthonormal set and 

ll/ll" = D(/(«).<’/.(')>P. (H-40) 

k 

which is the the Parseval’s theorem. 

There are different categories of wavelets, orthogonal wavelets, non-orthogonal, 
and biorthogonal wavelets. The Daubechies family, Symmlet, Coiflet, and Meyer 
wavelets are examples of orthogonal wavelets, while the Morlet wavelet is an example 
of a non-orthogonal wavelet. As mentioned earlier, not all applications of signal 
expansions require perfect reconstruction. Applications such as signal identification 
and classification are examples of this type. Therefore, unless the orthogonality is 
required to meet another criteria, we can exploit non-orthogonal wavelet transforms. 

2. Continuous Wavelet Transform 

The continuous wavelet transform (CWT) forms the mathematical basis for 
wavelet analysis. The concept behind wavelet analysis is that all basis functions can 
be generated from a single function called the mother wavelet, usually denoted by 
V’(t)- Other wavelets can be generated using two distinct operations; scaling and 
translation. Scaling means expanding (dilating) or compressing the wavelet function 
according to a specific scaling value. The scale will be denoted by s. The translation 
allows shifting of the (scaled) wavelet to a desired position in time. This shift will be 
denoted by k. The scaled and translated wavelet is denoted by 

Mt) = . (11.41) 

where y/s is a normalization factor to maintain the norm of '>ps,k{t) fo be constant 
at any scale s . 
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The integral form of CWT of the finite energy signal f{t) with respect to the 
wavelet function ip{t) is given by [Fief. 42, 46] 

/ OO 

fmsAt)dt, (11.42) 

-OO 

or 

W,(s, k) = ~ /_“ mi> *■ (11.43) 

The wavelet analysis is carried out by computing inner products (projections) 
between the signal and the wavelet functions. We can also interpret the wavelet anal- 
ysis as a linear operation which transforms the signal using modified kernel functions, 
where the kernel of the transform is the mother wavelet, and the modifications are the 
scaling and the translation operation [Ref. 44]. If the wavelet function satisfies the 
admissibility condition, we can apply the rule of the resolution of identity to recover 
the signal from its wavelet transform. The admissibility condition is given by [Ref. 
44] 

/•OO |ijr(a;)p 

Q = / ' \ /' dw < OO, (11.44) 

J — OO 

where ’J'(o;) is the Fourier transform of 'ipit). This condition will be true if ^(0) = 0 
or equivalently 

/ OO 

V’(t)dt = 0. (11.45) 

-OO 

The admissibility condition implies that the wavelet has a zero mean (i.e., DC com- 
ponent is zero). Therefore, it can be interpreted as a bandpass function. This implies 
that the wavelet must be an oscillatory function. The recovery (synthesis) formula, 
extracted from the resolution of identity rule, is given by 

rco roo 

f{t)= / Wf{s,k)'tl)s,k{t)dsdk. (11.46) 

Joo Joo 

3. Properties of Continuous Wavelet Transform 

The following properties of the wavelet ransform are given below without proof. 
More details can be found in [Ref. 45]. 
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1. Linearity: The wavelet transform is a linear operation on the analyzed signal. 
If the signal f{t) is given as: 

/W = AW + AW> W-47) 

then the wavelet transform of f{t) is given by 

Wf{s, k) = Wf,{s, k) + Wf,{s, k). (11.48) 

2. Shift-invariance: The CWT is shift invariant. If g{t) = f{t — to), then it has 
a CWT given by 

Wg{s,k)=Wf{s,k-to). (11.49) 

3. Time-scaling-invariance: If g(t) = ■^/(f )) if a CWT given by 

W,{s,k) = W,(^^A. (11.50) 

4. Parseval’s relation for the CWT becomes: 

r\\mfdt=^rr\ms,k)f^ nni) 

Joo Joo Joo S 

implying that the CWT conserves the signal energy. 

5. Time localization: Prom the definition of CWT and by the sifting property of 
the Dirac delta function, the CWT of an impulse at time to will be [Ref. 47] 

This means that the response of the CWT, at any scale s, will be 
a scaled and time reversed replica of the mother wavelet centered at location 
to on the shift axis (k). Therefore, the ability to define the location to will 
depend on the width of the scale wavelet. And since, a smaller scale value will 
result in a shorter wavelet, the time localization will be more precise. 

6. Frequency localization: In fact, the uncertainty principle controls the time 
and frequency resolution of the wavelet transform at any scale, as the smaller 
the scale is the shorter the wavelet and the wider its frequency representation 
are. Consequently, the shorter wavelet has worse frequency localization than 
the original wavelet. Generally, the time localization is better at smaller scale 
values while the frequency localization is better at larger scale values. 

4. Scalogram 

Using a scalogram, the finite energy signal f{t) is represented by the distribu- 
tion of \Wf{s, k)]"^ over the time-scale plane. From Parseval’s theorem for the wavelet 
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transform, Equation 11-51, integration of \Wf{s,k)\‘^ with respect to the differential 
provides the total energy. By the definition of the scaled wavelet function in 
Equation 11.41, the scale value s is inversely proportional to the frequency. Hence, 
the differential ^ is proportional to the differential of the frequency. Therefore, the 
quantity \Wf{s, A:)p may be viewed as a spectral density in units of power per Hertz 
[Ref. 44]. Consequently, the scalogram represents the power spectral density of the 
signal over the time-scale plane. Finally, we have the quantity 

(n.52) 

which represents the signal instantaneous power, while the quantity 

r m{s,k)W‘dk (11.53) 

Joo Joo 

represents the portion of the signal energy contained within the scale s. This fact 
is exploited in identifying the scale of each hop. In summary, the CWT is a linear, 
time-shift-invariant, time-scaling-invariant, and frequency-scaling-invariant operator. 

5. Discrete Wavelet Transform 

The CWT is defined by an integral transform over continuous variables in 
the scale s and the time shift k. Hence, performing the CWT requires expensive 
computations. Therefore, in practice a discrete grid for s and k is used. A widely 
accepted discretization is to specify s = and k = nkgs’^, where m and n are 
integers, and Sg > 1, and A:,, > 0 [Ref. 45]. Furthermore, if we select Sg = 2 and 
kg •= 1 we obtain the well-known dyadic wavelet sampling (tiling) grid. Hence the 
scaled and translated wavelet indexed by m and n is given by 

■0m, n (t) = 2^ 0 - n) . (II. 54) 

Figure 1 shows the time-frequency tiling of the STFT while Figure 2 shows the time- 
scale tiling of the wavelet transform. 
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In digital signal processing there are two related concepts that help under- 
standing the wavelet transforms: the multiresolution analysis and the filter-bank 
theory. These concepts are discussed next. 

6. Multi-Resolution Concepts in Discrete Wavelet 
Transforms 

Introduction. We will briefly introduce the WT from the perspective of the mul- 
tiresolution analysis (MR). The signal (or the time function) f{t) is expanded in 
terms of the wavelet functions. These wavelets have a frequency bandpass shape, so 
they result in a set of successive details of the signal. For the approximation we need 
a special basis, called the scaling function c^(t), which is not a wavelet. It has a low 
pass frequency behavior and performs averaging. The discrete wavelet transform is 
given by 

OO 00 OO 

/(t)= E c{k)Mt) + T. I: <i(j,k)i,d,k)(t), (11.55) 

k=—oo j=0/s=— OO 

where j and k are integers, the coefficient c{k) constitutes the coefficients of the 
approximation, while d{j, k) constitutes the coefficients of the added details or equiv- 
alently the fine resolutions [Ref. 47]. If the wavelets and the associated scaling 
functions form an orthonormal set of basis functions the coefficients are given by 



c(^) = 


(11.56) 


d{j, k) = {f{t),'(pj,k{t))- 


(11.57) 



The expansion form of the DWT is related to the subspace representation. Let 
us define the following subspaces: Vj is the scaling space at the level, and Wj is the 
wavelet space at the level. Let Vj+i = Vj+Wj. This means that the approximation 
at the {j 4-1)®^ level (or scale) can be represented by a coarser approximation Vj and 
a coarser detail Wj, both at the level. Note that the subspaces Vj are spanned by 
the scaling functions (f>k{i) and the subspaces Wj are spanned by the wavelet functions 
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Subspaces have to satisfy a set of conditions to become an MR representation. 
These conditions are [Ref. 46]: 

1. Orthogonality and completeness: 

VjDWj = { 0 }. 

00 

= L\ (11.58) 

3=0 

2. Scale-invariance: If f{t) £ Vj then f{2t) € Vj+i. 

3. Shift-invariance: If f{t) G Vj then f{t — ko) G Vj for a given constant ko . 

Scaling and Wavelet Equations. The MR subspace representation leads to a 
method to formulate two equations in terms of the unknown scaling and wavelet 
functions. Prom the orthogonality and completeness property we deduce that Vq C 
Vi. This means any function living in Vo, e.g., can be expressed as a linear 

combination of functions living in Vi, e.g., Hence, by substitution we get 

= V2Yi ho{k)(l){2t - k), (11.59) 

k 

where (poi’f-) = and (f>i{t) = <f>{2t). This equation includes two different scales of 
the scaling function and is known by many names: the scaling equation, the dilation 
equation, the refinement equation, or the multiresolution-analysis equation. The 
coefficients ho{n) are called the scaling filter coefficients. Using the condition, Vq © 
Wo = Vi, which implies that Wq C V, we have 

^(t) = ^/2Y,hl{k)(f>{2t - k) (11.60) 

k 

also is called the wavelet equation. The coefficients h\{n) are called the wavelet filter 
coefficients. 

In summary, an MR approach leads to a complete orthonormal set of the 
scaling and wavelet functions. Moreover, if the scaling function has compact support. 
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the wavelet function will be built up of a finite number of the scaling multipliers. 
Therefore, the scaling and the wavelet functions can be realized by finite-impulse- 
response (FIR) filters. For more details about wavelet and scaling equations see [Ref. 
47, 42]. 

For discrete data the filter-bank concept leads to a simple method for com- 
puting the wavelet coefficients. The wavelet function is replaced by the coefficients 
of the wavelet filter hi{n) and the scaling function is replaced by the coefficients of 
the scaling filter ho(n). In [Ref. 46, 47, 42], the following two recursive equations are 
obtained: 

k) = - 2k)c{j + l,m) 

m 

d{j, k) = X! hi(m - 2k)c{j + 1, m). (11.61) 

m 

These two recursive equations enable us to compute the j scale wavelet trans- 
form from the (j + 1)** scale using ho{m) and hi{m). The factor 2 in front of the 
parameter k decimates the output wavelet coefficients by 2. Figure 3 shows the re- 
alization of the recursive equations. This realization is equivalent to the well-known 
2-band analysis bank of filter implementation. The implementation of the discrete 
wavelet transform using the recursive procedure is known as Mallat’s algorithm. Fig- 
ure 4 shows a multi-stage wavelet analysis bank and the corresponding ideal subspace 
designation (or spectral decomposition). 

7. Daubechies Wavelet Family 

The Daubechies wavelet family is a compactly supported orthonormal set of 
wavelet functions [Ref. 42]. The Daubechies wavelet are generated by solving the 
scaling and wavelet equations. An additional set of constraints is applied to satisfy the 
maximum number of vanishing moments for each wavelet. Each Daubechies wavelet 
is assigned a number related to the number of vanishing moments, hence, each wavelet 
is denoted as Daub-N. There are two conventions in assigning the number N; the first 
one is to let N refer to the number of vanishing moments; the second one is to let N 
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refer to the length of the wavelet filter (i.e., the number of wavelet coefficients). We 
will adopt the first designation since it is adopted by the MATLAB Wavelet Toolbox, 
where each Daubechies wavelet is expressed by the coeflScients of its scaling filter. 
The coefficients of the wavelet filters are the interpolated high pass versions of the 
scaling filter. The Daubechies wavelet of order N has 2N coefficients, and has finite 
support over [0, 2A^ — 1]. The number of vanishing moments is an indication of the 
smoothness of the the wavelet filter, since the number of vanishing moments will 
imply the number of zeros of at w = tt. The higher the order the longer and 
smoother is the Daubechies filter. 
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Figure 1. Time-Frequency tiling of the STFT. 
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Figure 2. Time-Scale tiling of the wavelet transform. 
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Figure 3. Recursive computation of wavelet coefficients. 
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Figure 4. Multiple-Stage wavelet analysis bank and its subspace designation. 
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III. WAVELET TRANSFORMS AND 
CORRELATION FUNCTIONS 

Correlating two functions results in a measure of similarity between them 
and is accomplished by evaluating an inner product. Correlation analysis has been 
used extensively in the signal processing and communication area, e.g., in spectral 
estimation, system modeling, signal detection and parameter estimation. 

The Wiener-Khinchin theorem relates the signal’s auto-correlation function 
and power spectral density for a stationary stochastic process, through a Fourier 
transform. Due to their time- varying spectra, FH signals are non-stationary and can 
be represented using time-frequency distributions. Traditionally, Fourier kernels are 
used in the time-frequency distributions. 

Wavelet decomposition can be used to represent non-stationary signals over the 
time-scale plane. Therefore, we will consider the wavelet transform of the correlation 
function as an alternative for non-stationary signal representation. 

In this chapter, we will introduce different choices for the correlation function 
and derive the corresponding wavelet transform response. In the derivation of the 
wavelet transform we will use the analytic definition of the correlation function and 
the wavelet function without specifying a functional form. 

A. CORRELATION FUNCTIONS 

1. Definitions of Auto-Correlation Functions 

Depending on the underlying process, various definitions are given to the 
auto-correlation function (ACF). The process may be deterministic, stochastic, finite- 
energy, infinite-energy, non-time- varying (stationary) or time- varying (non-stationary). 
Next, we will introduce the different definitions of the ACF. The main reference for 
this topic is [Ref. 51]. 

1. Deterministic ACF: The ACF is defined for deterministic, infinite-energy, and 
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non-time- varying signals by 

1 T 

R{t') = ^ + r)dt , (III.l) 

where "T" is the observation interval and "r" is the time lag or the time delay. 
For finite energy signals the ACF is defined as: 

/ OO 

x{t)x* {t + r)dt . 

•OO 

For stationary signals or processes (i.e., wide-sense stationary) R{t,r) is a 
function of the time delay only, hence, it is denoted by R{t). For a non- 
stationary signal, R{t, r) will be represented by a 2-D surface over the time 
"t" and the time delay "r". 

2. Stochastic ACF: The ACF of a stochastic process is the correlation between 
two samples of the process taken at ti and ^ 2 , and is defined as: 

= (in.3) 

where E is the expectation operator and stands for the complex conju- 
gation. For a stationary (or wide-sense stationary) process, R{ti,t 2 ) depends 
only on the time lag r = resulting in a stationary spectrum. For a non- 

stationary process, R{ti,t 2 ) depends on the time instants and ^ 2 , resulting 
in a time- varying spectrum. 

2. Properties of Auto-Correlation Functions 
We will briefly introduce the significant properties of the ACF for finite-energy 
signals. For more details see [Ref. 36, 52]. 

1. Conjugate symmetry: 

R{t,T) = R{t,-ry, . (III.4) 

i.e., the real part of i?(t, r) is an even function with respect to the lag r while 
the imaginary part is an odd function. 

2. Mean-squared value: 

= (in.5) 

3. Positive definiteness: The ACF, denoted by R{ti,t 2 ), is said to be positive 
definite if 

^aia*R{ti,tj)>0 (III.6) 

for any ai, ay ^ 0. Positive definite property of the ACF guarantees that the 
spectrum of the signal is non-negative. 
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4. Wiener-Khinchin theorem: 

/ OO 
“OO 

The power spectral density function Sxx{f) of the stationary signal is ob- 
tained by Fourier transforming its (positive definite) ACF. Note that the non- 
stationary signals have time-varying spectra, therefore, no power spectral den- 
sity is defined. 

3. The Instantaneous Correlation Function 

The ACF of deterministic and stochastic processes are computed using time 
domain averaging and the expectation operator, respectively. This means that a 
smoothing process has to be applied to compute the correlation functions. 

The instantaneous correlation function (ICF) does not use an averaging oper- 
ation (integration nor expectation). The instantaneous correlation function is defined 
as the product of two samples of the process. These two samples are drawn at two 
time instants centered about time t. The instantaneous correlation function R*{t,r) 
is defined as: 

R\t, r) = a: + 0 x* (’f “ 0 » (HI.8) 

where "i" stands for the instantaneous nature of the correlation function. 

R'{t,r) satisfies the following properties: 

1. Conjugate symmetry, i.e., 

R\t,T) = R'\t,-T). (III.9) 

2. Squared value, 

R'■{^,T)|.=„=|a:{i)|^ (III.10) 

i.e., at r = 0 , R^{t,r) will represent the instantaneous power of the signal at 
time “t”. 

If x{t) is a sinusoidal signal then multiplication of the instantaneous values of 
x[t) in the sense of Equation III-8 will generate cross terms in the ICF. For example, 
the real-valued sinusoidal signal x{t) = A cos(wt) has an ACF given by: 

R{t) = -^cos(wr), (III.ll) 
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while the ICF is given by: 



R*{t, t) = — [cos(2c<;t) + cos(o;t)] . (III. 12) 

The ACF of a single sinusoidal signal has only one component and no cross 
term, while the ICF has cross terms. If the signal x{t) is represented by its analytic 
form, say x{t) = Aexp jwt, then its ICF is given by: 

R*{t,r) = exp {jur). (III.13) 

That is, the ICF of a single complex exponential signal has no cross term. Therefore, 
to minimize cross terms it is recommended to use the analytic form of the input 
process in the computation of the Time- Frequency Distributions [Ref. 49]. In this 
dissertation, we will focus on the analytic form of the signals. Therefore, the first 
step of the processing scheme will be transforming the intercepted (real) signal using 
the Hilbert transform technique. 

Another point worth mentioning is the exploitation of the conjugate-symmetry 
property of the ICF. In computing the surface of the ICF we need to compute only 
half the ICF surface, corresponding to positive (or negative) time lag “r". The other 
half can be generated by the conjugate-symmetry property. This implies that the 
second half of the surface of the ICF carries no additional information. Therefore, 
the processing scheme may be applied to only one half of the ICF surface which will 
save computational cost and storage. 

B. WAVELET TRANSFORMS OF CORRELATION 
FUNCTIONS 

In Chapter II, we discussed the linearity property of the wavelet transform 
which allows the application of linear system theory to the wavelet transform. In this 
section, we will investigate the wavelet transform for the different definitions of the 
correlation function. 
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1. The Wavelet Transform of the Auto-Correlation 
Function 

The wavelet transform of the stochastic ACF of a stationary process is ad- 
dressed in [Ref. 48], similarly, we will address the ACF for a deterministic signal. 
Wavelet transform of the ACF will have similar results. For a deterministic, finite- 
energy, stationary signal, the (time-averaged) ACF is given by: 

/ OO 

x{t)x* {t + T)dt . (III. 14) 

-OO 

From the Wiener-Kh inchin theorem we get 

/ OO 

R(r)e-^^-f-^dr. (III.15) 

-OO 

Let Wxx{s,t) denote the Wavelet transform of R{t). Note, the subscript “rex” 
in Wxx{s,i) stands for the wavelet transform of the ACF of x{t) (in contrast to 
Wx{Sil) which denotes the wavelet transform of a;(<)). The wavelet basis function 
is denoted by g{r). The wavelet transformation will be carried over the time lag 
variable “r” to the wavelet shift variable “t and the wavelet scale “s”. Wxx{s^i) is 
then given by: 

Wxx{s, = RW 

Let G{f) denote the FT of p(r), then 

FT{g (^)} = kl . (IIU7) 

For positive scale values , “s” , g is given by: 

= r . (IIL18) 

S J J—oo 

We can write Wxx{s,i) as: 

rOO roo 

Wxx{s,e) = V^ / R{r)e-^^^^^drG*{sf)e^'^^^^df, (III.19) 

J — OO •/—OO 

or 

roo 

Wxxis, i) = >A Sxx{f)G*{sf)eP^^Rdf . (III.20) 

J—OO 
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Equation III. 20 has the form of an inverse FT from the variable / back to the 
variable £. Therefore, we can write 

w„{s, e) = FT-‘{ (Vis..(/)G'(S/))} . (III.21) 

From Equation III. 21 we deduce the following: 

1. The wavelet transform, of a finite-energy signal, at any wavelet scale s ^ 0, 
represents a linear filtering operation using a band pass filter whose impulse 
response is the (time-reversed) wavelet function at the scale s. Equivalently, 
the filter has a frequency response given by the FT of the scaled wavelet. 

2. The wavelet transform of the ACF, R{t), of the stationary finite-energy signal 
x{t), gives a band pass filtered version of the power spectral density Sxx{f) of 
this signal (up to a constant, y/s), the used band pass filter is dependent on 
both the chosen wavelet function and the chosen wavelet scale. 

2. The Wavelet Transform of the Instantaneous 
Correlation Function 

The Wigner-Ville Distribution (WVD) is used to represent non-stationary pro- 
cesses since they are characterized by their time- varying spectra. The WVD applies a 
one-dimensional Fourier transformation to the ICF. The Fourier transform is carried 
out taking the time delay r to the frequency /, keeping the global time t unchanged. 
This allows for displaying the time evolution of the spectrum of- the signal. For 
one-dimensional time signals, the one-dimensional wavelet transform carries out a 
transformation from one global time variable “t” to the two wavelet variables, the 
shift £ and the scale s. Consequently, the signal is represented by a time-scale distri- 
bution in the wavelet domain. The wavelet domain is called the time-scale domain. 
For the two-dimensional surface, indexed by time “f ’ and the delay “r” , we carry out 
the wavelet transformation along the time delay. This permits a display as a func- 
tion of time. In this section, we will derive the relation between the time-frequency 
representation of the signal x{t) using the WVD and the time-scale representation of 
the same signal using the wavelet transform of the ICF. 
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Let Vx{t, f) denote the WVD of the signal x(t), 

Since the ICF of the signal x{t) is defined as: 

(III.23) 

where r is time shift (or time delay) relative to time t. The WVD may be defined as: 



/ OO • , 

Rl{t,r)e-^^^f^dT. (III.24) 

-OO 

The wavelet function g{r), has been scaled by s and shifted by i along the r axis, to 
form gs,i{T) which is given by: 

g,Ar) = 

Let G{f) denote the FT of g{r), thus 

FT{g } = s G{sf)e -^^^^^ . (III.26) 

Then 

= r sG{sf)e-^^^^^e^^^^^df . (III.27) 

S J J—oo 

Let W*j.{t;s,£) be the wavelet transform of the ICF R].{t,T). Note that the 
superscript ‘T’ indicates the instantaneous feature of the ICF. s, i) is given by: 

^ WUt; s. Q = :;)j /_“ r)g' dr . (III.28) 

Substituting g* from Equation III.27 into Equation III.28 and exchang- 
ing the integration operations we get 




for s > 0 . 



(III.25) 



W4(t,s,£) = x/i r G*{sf) [ r Ri{t,T)e-^^-f^dr 

J—OO U—oo 






Substituting Equation III.24 into Equation III. 29 we have 

roo 

WUt; sJ) = V^ G'{sf)V,{t, . 

J—OO 



(111.29) 



(III.30) 
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Equation III. 30 is in the form of an inverse Fourier transform. Therefore, s, i) 

and \^sG* {sf)Vx{t, f) are a Fourier transform pair with respect to the variable 
i.e., This relation suggests that we can obtain a filtered version of the WVD by 
Fourier transforming the W^xit; s, £). 
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IV. FREQUENCY HOPPED SIGNALS AND 
THEIR CORRELATION FUNCTIONS 



Communication systems can utilize a large number of digital modulation tech- 
niques. Among those, spread spectrum modulation is widely used. Spread spectrum 
refers to any modulation scheme that produces a transmitted bandwidth much larger 
than the information bandwidth [Ref. 50]. We will briefly address the different digital 
modulation schemes and focus on frequency hopping. 

A. DIGITAL MODULATION SCHEMES 

For comparison, digital communication schemes are briefly presented. Digital 
modulation techniques use the binary and the M-ary schemes. Binary schemes consist 
of the on-off keying (OOK), also called ASK, binary phase shift keying (BPSK), and 
frequency shift keying (FSK). M-ary schemes are generalizations of the binary schemes 
for transmitting M symbols, e.g., M-ary PSK or M-ary FSK. 

1. Binary schemes 

The major reference of this topic is [Ref. 51]. 

• OOK: The OOK signal is represented by: 

s{t) = Ac m{t) cos Wct^ 

where Ac is the carrier signal amplitude, m{t) is the transmitted binary data 
over the bit duration Ti,. For unipolar representation it has either “1” or “0”. 
Therefore, over an observation time Tj the signal is zero (i.e., m{t) = 0) or it 
is a single sinusoid at frequency Wc (i.e., m{t) = 1). 

• BPSK: The BPSK signal is represented by: 

s{t) = Ac cos [wct + A^m(t)] , (rV.2) 

where m{t) is the bi-polar message signal. Over the bit interval Tj, the signal 
has the value ±1 depending on whether the bit information is 1 or 0. A9 is 
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the corresponding phase shift. Usually, A9 is chosen as tt/2. Therefore, the 
BPSK signal has either one of the two forms: 



s{t) = dtAcSinWct ; for m(t) = ±1 . 



(IV.3) 



Therefore, under any value of m{t) the BPSK signal will be a sinusoid with a 
fixed frequency component but with a phase of zero or tt. 

• FSK: The continuous phase FSK signal is given by: 



where Df is a, modulation index. For binary messages m{t) the resultant 
FSK signal is called the binary FSK scheme. The FSK signal has one of two 
frequencies without phase discontinuities at the transition points. 

2. M-ary Schemes 

For an M-ary schemes a message m{t) has M symbols. Consequently, the 
transmitted signal will have M different states. 

• M-ary ASK: The M-ary scheme of the OOK (or ASK) may be implemented 
for different values of M. An example is QAM (i.e., M = 4). All of the QAM 
states have the same single frequency components but differ in amplitude. 

• M-ary PSK (MPSK): The M-ary PSK is generated similar to BPSK with the 
exception that the value of Ad is chosen according to the number M. Thus, 
all the MPSK states have the same frequency component but differ in phase. 

• M-ary FSK (MFSK): The M-ary scheme is similar to BFSK. But, instead of 
two symbols (states), it has a set of M symbols. Consequently, M different 
frequencies are transmitted. Thus, the MFSK signal can be expressed as s(t) = 
AcCos{wit + 6i) where Ui and 6i are the frequency and phase corresponding 
to the symbol of the message m{t). 

In summary, digital modulation techniques are characterized by their signaling 
mode which consists of a single fixed frequency component over the duration of a given 
bit. 



s{t) = Ac cos [Wct -t- 0(t)] , 



(IV.4) 



and 




(IV.5) 
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B. SPREAD SPECTRUM COMMUNICATION SIGNALS 

The major reference for this section is [Ref. 19, 50]. Spread spectrum (SS) 
communication signals are characterized by a wide transmission bandwidth, and a 
low power spectral density. SS signals have two main advantages; 

1) The message has a low probability of being intercepted (LPI) as a result of 
the wide frequency band, and the low power spectral density of the signals. 

2) SS systems can reject jamming signals and allow many users to share the same 
frequency band as a result of the spreading gain [Ref. 19, 50]. 

Among the different possible SS modulation formats, the follovdng are preva- 
lent: 

1) Frequency Hopping (FH): The complex baseband signal, c{t), with basic pulse 
shape p{t), is given by 

c(i) = [i (27r/nt + (l>n] p[t - nTh], (IV.6) 

where Th is the pulse duration, better known as the hop interval. The pseudo- 
randomly generated sequence of frequencies {/„} will drive the modulator to 
generate a modulated version of p{t). {4>n} is an associated phase shift at each 
carrier frequency /„. 

2) Direct Sequence (DS) Modulation: The complex baseband signal, c(t), of DS 
is given by 

c{t) = (IV.7) 

n 

where {cn} is a pseudorandom sequence which modulates a sequence of pulses 
over a duration Tc, known as the chip interval. 

3) Time Hopping (TH): The pulse waveform is given a fraction of duration Ts, 
i.e., Tg/M. A typical time hopping waveform might be 

cW = Ep(*-(" + ^)r.). (IV.8) 

which means the pulse will be controlled by the pseudo-random number (o„) 
to appear at any of the M time segments within the duration T^. 

4) Hybrid Modulations: A blend of the above techniques may satisfy better per- 
formance according to some requirements. 
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Figure 5 shows a typical SS system. The SS spectrum uses two types of mod- 
ulation. The first type, the baseband data modulation, is also called the primary 
modulation. The second type, the actual SS modulation, is also called the secondary 
modulation. For simple SS system realizations, certain combinations of primary and 
secondary modulations are usually employed. For DS, the binary phase shift keying 
(BPSK) is used as the primary data modulation scheme. For FH, the M-ary frequency 
shift keying (MFSK) is used. The FH system which implements MFSK data modula- 
tion is known as the pure FH scheme, and is the most popular and widely applied FH 
scheme. Another advantage for the pure FH scheme is its resistance against Repeat- 
back jammers (RBJ) which estimate the FH signal frequency and consequently send 
a proper jamming waveform. Thus, pure FH schemes minimize the hostile activity of 
RBJ as each hop frequency carries one symbol of the transmitted message only. 

C. THE FREQUENCY HOPPED SIGNAL 

The FH of the BFSK signal is given by 

x{t) = y/‘lP sin [(wq "b "b ■> -9) 

where P is the signal power. The frequencies uq and Ao; are constants, and is the 
frequency for the symbol dn whose value is ±1. 

During any hop, only one frequency component will be at the transmitter 
output when the hopping interval equals the symbol interval. The range of variation 
of ujc = iojn -b o;o ± Aw) is known as the hopping bandwidth. Two types of FH 
systems exist depending on the hop interval and the symbol interval T^: fast- 
frequency hopping (FFH) and slow frequency hopping (SFH) which differ in their 
hopping speed. The number of hops within one symbol duration is the measure of 
the speed of the hopping rate. For pure FH signals with BFSK, the SFH requires 
Th > Tg while FFH requires < Tg. The SFH contains one or more symbols during 
the hop interval while the FFH has one or more hops over one symbol interval. Figure 
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6 shows a typical pure FH system using the BFSK as a primary modulation. Figure 

7 illustrates a typical time-frequency behavior of the BFSK pure FH scheme. 

D. THE INSTANTANEOUS CORRELATION FUNCTION 
OF FREQUENCY HOPPED SIGNALS 

Spread spectrum studies have usually considered the FH signal as a stationary 
process [Ref. 19, 50] even though the spectrum of the FH signal varies with each 
hop interval. Therefore, the correlation representation, using time averaging, is not 
suitable for this nonstationary process. 

One way to identify the FH signal is to monitor the time evolution of the signal. 
Hence, we need to keep the time dependency in the correlation representation. In 
this work we resort to the time-varying correlation definition for characterizing the 
FH signals. Therefore, we select the instantaneous correlation function (IGF) as the 
candidate correlation representation. In this section, we will address the structure of 
the IGF of the pure FH signals. 

The pure FH signal may be represented as successive intervals (i.e., hops) 
with single complex exponentials. The frequency within each interval (i.e., the hop 
frequency) is controlled by a randomly generated sequence of numbers, therefore, we 
can assume the following about the used frequencies: 

(1) The number of the different hopping frequencies is much larger than the num- 
ber of observed hops. 

(2) The selection of discrete hopping frequencies (hopping pattern) is determined 
by a pseudo-random number generator (PNG), i.e., all hopping frequencies are 
equi-probable. 

(3) The observation interval is much smaller than the period of PNG. The period 
of PNG is the number of generated hops times the hop interval. 

As a result, we can deduce that any two successive hops will have two dif- 
ferent frequencies. The difference in the frequency of adjacent hops will generate 
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the structure of the instantaneous correlation functions. Recall that we define the 
instantaneous correlation function as: 

= + . (IV.IO) 

where t is time and r is time lag. Over the two-dimensional plane of time and time 
lag (i.e., t — T plane) we compute the IGF for the time lag |r| < Th. which will allow 
correlating adjacent hops only. If the values of {t -f |) and {t — |) are both confined 
within the hop then the IGF will be given by: 

R{t,r) = 

= , (IV.ll) 

which is a function of r and ul only. Note that the values of {t -I- and {t — |) are 
both confined within the same hop if they satisfy, 

(L - l)Th < (< + ^) m. (IV.12) 

This inequality forms the boundaries of the diamond cellular shape for different 
values of L. This cellular structure is- shown in Figure 8. Inside each diamond Di, 
the IGF results from correlating signals belonging to the same hop, while outside the 
diamond the IGF results from correlating signals belonging to two consecutive hops. 

A detailed derivation of the IGF is given in Appendix A, which leads to the 
doubly indexed correlation function Rm,a{t,T) given by: 

'^) = exp j I (Um ~ -|- (u>m + • (IV.13) 

where m and n are the indices of the two adjacent hops. Note that: 

1) Within the main diamond of the n*** hop, i.e., m = n, the IGF is given by 

= exp{j(o;„r)}. (IV.14) 
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2) Outside the main diamond in the upper triangle between hop numbers (n) and 
(n + 1), the IGF is given by 

r) = exp + {un+i + ^n)^| • (IV.15) 

In conclusion, the instantaneous correlation function of a pure FH signal has 
a cellular or diamond structure. Inside the diamond it has a single complex expo- 
nential component along the delay axis representing the hop frequency. Outside 
the diamond, R{t,r) is a product of two terms, and where 

and Um are two consecutive hop frequencies. 
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Figure 5. Modern SS system. 



46 




Figure 6. BFSK-FH system. 
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Figure 7. Time-Frequency behavior of BFSK-FH signal. 
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Figure 8. FH signal and the cellular (diamond) structure of its ICF. 
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Figure 9. Areas and indices of the doubly indexed function Rm,n{t,T). 
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V. 



ANALYSIS USING MORLET WAVELETS 



We studied the structure of the ICF surface obtained for complex FH signal 
in chapter IV, and showed that it consists of complex sinusoidal components. In this 
chapter, we will analyze the Morlet wavelet transform of the ICF. We will use the 
Morlet basis function for two reasons. First, the Morlet wavelet is a complex sinusoid, 
modulated by a Gaussian window, and is inherently best suited for filtering sinusoidal 
signals. Second, the Morlet wavelet has a mathematical formulation which makes the 
analysis and derivations tractable. 

Recall that the mother Morlet wavelet g{t) is given by: 

g{t) = for — oo <t < +oo , (V.l) 

where k is a constant that represents the modulation frequency, and cr^ is inversely 
proportional to the roll-off factor of the Gaussian window. The Fourier transform is 
given by: 

G{u) (V.2) 

The Morlet wavelet does not satisfy the admissibility nor the orthogonality 
condition as it exists over an infinite time interval, but for practical applications it is 
truncated when it is sufficiently attenuated. If the wavelet has finite support (non-zero 
only over finite interval), then the scaled and shifted wavelet, gs,e{t), will be non-zero 
over the interval [i — nas, £ + nas], where s is the scale and i is the shift. Here, n 
is ai preselected number which ensures sufficient decay. The functional dependence of 
9s,e{t) is as described in Equation 11.41. Usually, n is chosen > 4 [Ref. 46]. 

A. TRANSFORM OF A TIME-LIMITED COMPLEX 
EXPONENTIAL 

We showed in Chapter IV that the ICF of the FH signal consists of complex 
exponential components with frequencies that differ from region to region. Therefore, 
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we need to address the response of the wavelet transform to complex exponential 
signals. Let the signal f{t) be a complex exponential defined over the interval 

t 6 [LT, (L + 1)T], where T is the hop interval, and L stands for the hop. When 
(7 = 1, the Morlet wavelet transform (MWT) is given by: 




£-^ns 

—ns 






(V.3) 



where £ is confined to the interval [LT + ns, {L + 1)T — ns]. This interval guarantees 
that the wavelet will be confined over the same hop without crossing the border to 
another hop. Later on we will compute the wavelet transform when the signal runs 
over the border between two adjacent hops. The exponent in Equation V.3 can be 
written as: 



/ _k e\ M 



Substituting Ul — k/s = u we obtain 



(V.4) 



W(s,e) = f‘*^exp 

V S J£—ns 



2s2 






dt . (V.5) 



This equation can be reduced to 



W(s,^) = A^yexp 



<j2,,2 jl2 



(EF), 



(V.6) 



where EF is given by: 



EF = erf 



1 



^{n + j{iVLS-k)) 



+ erf 



■^(n- - *:)) 



(V.7) 



for £ G [LT -|- sn, {L 1)T — sn]. 

We note that the magnitude of W (s, £) is independent of the wavelet shift 
variable £. 

The term EF defined above can be approximated using the following formula for the 
complex error function [Ref. 54] with 
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erf{x + jy) = erf{x) + - — [1 — cos 2a:y + j sin 2a:j/] 

^7i JU 



2 _ 2 « g-(l/4)m2 

+-e Z 2 ^2 , 4^2 Ifrnix, v) + JQraix, y)\ 

m=l ' ' 



+s{x,y), 



(V.8) 



where 



/m(2;, y) = 2a: — 2a: cosh (my) cos(2a:y) + m sinh(my) sin(2xy) 

g^(x,y) = 2xcosh(my) sin(2xy) + msinh(my) cos(2xy) , 

and |e(x,y)| « 10"^®|er/(x, y)|. 

Since EF is of the form 

EF{x, y) = erf{x + jy) + erf{x - jy), 

we need to determine erf{x — jy). Recall that sin(x) and sinh(x) are odd functions 
while cos(x) and cosh(x) are even functions. Thus we can express 



erf{x - jy) = erf(x) + [(1 - cos 2xy) - j sin 2xy] 



(V.9) 



2 _^2 ^ 

m^ + 4x^ 

m—1 



ifm{x, -y) + j9m{^, ~y )) , 



where 



fm{x, -y) 

9m{x, ~y) 
Or, equivalently, 



2x — 2x cosh(my) cos(2xy) + msinh(my) sin(2xy), 
— 2x cosh(my) sin(2xy) — msinh(my) cos(2xy) . 

fmiFt y) — fmip'iy)i 
9m{x,-y) = -9m{x,y). 
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Using 



erf{x -jy)= erf*{x + jy), 



(V.IO) 



and 

erf{x + jy) + erf{x - jy) = 2Re {erf{x + jy)} , (V.ll) 



where Re{ } denotes the real part, we can express W{s,£) as: 



W (s, €) = Ay/2irs exp 






2 - Y + + jojL^ 



■Re I erf f ^[n 4- j(wL5 - k)] 



(V.12) 



for I € [Z/T" 4” ns, {L 4~ 1)X' — ns]. 

Equation V.12 shows that the transform has a linear phase response, $(o;), 
independent of the scale s, which is given by: 



$(o;) = uii. 



(V.13) 



Note that if the complex exponential signal has an initial phase shift, Bi, the phase 
of the output of the transform will change to 

^{u)) = {udABi). (V.14) 

The amplitude may be expressed as AB{s,u), where B{s,oS) is the wavelet 
gain given by: 

erf - *=)] 

(V.15) 

The gain of the Morlet wavelet transform B{s, oS) is plotted in Figure 10 as a function 
of scale and frequency. 

The spectral response of the wavelet is plotted for values of s = 0.5,1, and 
2. Over this dyadic grid the spectral response contains some regions of low response 
due to the narrow bandwidth of the Morlet wavelet filters. In many applications, 
discrete-scale wavelets are preferred. Note that, the scale grid needs to be sampled 



B{s, of) — V2 tts exp 



2, ,2 j ^2 

;r 4- suik 






Re 
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linearly (i.e., 5 = 1, 2, 3,4, . . . ) or the roll-ofF factor (a^) of Equation V.l has to be 
decreased when applying Morlet wavelets for full spectral coverage. 

Complex error function. The error function used in W{s,i) has a complex ar- 
gument and we use a computational method given in [Ref. 54, 56] to compute it. 
The error function is defined as 

erf{z) = 4 = re~^"du, (V.16) 

\/7r ^0 

while the complementary error function is defined as: 

erfc{z) = 1 - erf{z). (V.17) 

Let w{z) be defined as: 

w(z) = e~^^erfc(—jz). (V.18) 

Then, for a complex argument z, the error function is expressed as: 

erf{z) = 1 — e~^^w{jz). (V.19) 

The function w{z) is tabulated for some values of 2 : in [Ref. 54] and may be 
computed using the algorithm described in [Ref. 56]. 

B. ANALYSIS OF THE TRANSITION REGION 

The transition region is defined as that where the FH signal f{t) changes from 
one to another hop at time tr- Thus, the wavelet transform involves two different 
signal frequency components in the transition region. The transition region between 
interval L and interval L -1- 1 will be defined as ^ € [{L + 1)T — ns, {L + l)T + ns], 
where the Morlet wavelet support interval is 2ns. 

Hence, the Morlet wavelet transform is given by: 

= 71 L (¥) + Ts C (¥) 

where tr — ns<i<tr+ ns. 
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The two parts of Equation V.20 can be solved similarly to Equation V.3, which 



leads to 



w{s,e) = 



Pks 

/-exp 



s^(jo\ k'^ 7 • /) 

— ^ + su)i,k + 



erf 






tr-i 



Pks 

+VT*^ 



,2 



-j{<^LS-k) 



+ erf + j - A:)] 



S LOt,-, k , . „ 

2 2 ^<^L+ik + Jojl+i^ 



|er/ - j{^L+iS - k)^ 



-erf 



'tr - i. 



.V2 



- j{u)L+\S - k) 



(V.21) 



Over the length of the transition region the magnitude of MWT depends on 

t 

Figures 11 and 12 show an example of the magnitude and the phase of the wavelet 
transform over a transition region. Two frequencies are used (4 and 7 rad/sec) and 
the magnitude and the phase of the wavelet transform are plotted for scale s = 1. 
The magnitude of the wavelet coefficients, for the signal with the first frequency (i.e., 
f=4), is larger than the magnitude of the wavelet coefficients of the signal with the 
second frequency (i.e., f=7). The magnitude and the phase of the wavelet transform 
for scale s = 0.5 is plotted in Figure 12. The magnitude of the wavelet coefficients, 
for the signal with the first frequency (i.e., f=4), is smaller than the magnitude of 
the wavelet coefficients of the signal with the second frequency (i.e., f=7). Note also 
that there is a gradual transition between the two magnitudes according to Equation 
V.21. This indicates that the Morlet wavelet transform does not produce spikes at 
the transition points between frequency hops. 
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C. ANALYSIS OF IDEAL WHITE GAUSSIAN NOISE 



We investigate the response of the Morlet wavelet to ideal white Gaussian 
noise. The WT is given by: 

W„(sJ)= r -^n{t)g-{—)dt. (V.22) 

J—oo Y S S 

For white Gaussian noise n{t) and due to the linearity of the wavelet transform, the 
output Wn{s,i) has a Gaussian distribution. Since n{t) is assumed having zero mean 
and variance then the mean value of Wn{s,i) is 



E{Wn{s,e)} = E 



= 0 . 



(V.23) 



Since Wn{s,i) is zero mean, its variance is given by: 
d = E{\W{s,l)?} = E{W(s,e)W’{s,l)} 

= \ C II (^) (^) ) ■ 

(V.24) 



When two samples n{ti) and n{t 2 ) are assumed independent identically distributed 
Gaussian (i.i.d.) random variables, the term E{[n{ti)n* {t 2 )} is equal to zero except 
when ti = Hence we have 




(V.25) 



Consequently, the Morlet wavelet transform of white Gaussian noise, with zero 
mean and variance, is Gaussian with zero mean and a variance of 
wavelet scale s. Note that we assumed an infinite support for the Morlet wavelet in the 
above derivation. If we carry out the same procedure for the finite support wavelet. 
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we will have a different expression for the variance of Wn{s, £) due to integrating the 
term \g over a finite interval. In such a case, the variance becomes 



S J-ns I \ ^ / 



dt , 



which results in 



= al^/^erf{ns) 



(V.26) 



(V.27) 



We note that for ns » 1, erf {ns) is well approximated by 1. Thus, we conclude that 
the variance of the wavelet transform will be essentially independent of the wavelet 
scale and can be approximated by o\ypK. 

Actual Noise of the Wavelet Surfaces. Application of the wavelet transform 
to the ICF surface will result in a number of wavelet surfaces corresponding to the 
number of wavelet scales used. Since, we exploit the wavelet surfaces for identify- 
ing the FH signal, we need to investigate the noise distribution over these surfaces. 
The noise of the wavelet surfaces is approximated as Gaussian noise. The following 
considerations assist in making this decision: 



1) The noise background is assumed to be additive white Gaussian noise. The FH 
sinusoidal components and the noise are assumed to be independent. Noise 
realizations at two different time instants are assumed to be uncorrelated and 
independent. 

2) The ICF surface consists of noise components which result from the product of 
two independent Gaussian random variables. This resulting noise surface has a 

-distribution shape shown in Figure 13, where r is the correlation coefficient 
of the two Gaussian random variables [Ref. 57]. 

3) The wavelet transform of the ICF surface is a weighted sum of the ICF samples. 

According to the central limit theory, usually, the sum of number of identically 
distributed random variables tends to a Normal (Gaussian) distribution. Since we 
take the wavelet transform of the ICF surface in the time delay direction, the samples 
of the wavelet surfaces, in the time delay direction, may be approximated as having 
a Gaussian distribution. Simulation results show that the noise distribution over 
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the wavelet surfaces, taken in direction of time delay, has approximately a Gaussian 
distribution. The Gaussian assumption was tested using a chi-square test. The main 
central part of the noise distribution (excluding the tails) was found to be Gaussian 
with confidence level > 90%. 
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Morlet wavelet gain 




Figure 10. Morlet wavelet gain. 




angle( W(s,l) ) _ ■ ^V(s.l) ' 



s= 1 , w1=4, w2= 7 




time shift, I 




Figure 11. Morlet wavelet gain and phase over the transition region (s=l). 
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angle( W(s,l) ) W(s,l) I 



s= 0.5 , w1=4, w2= 7 





Figure 12. Morlet wavelet gain and phase over the transition region (s=0.5). 
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Figure 13. The K-distribution. 
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VI. 



ANALYSIS OF THE PROCESSING 
SCHEME 



A. INTRODUCTION 

Our goal is to exploit the wavelet transform of the ICF surface in order to 
identify FH modulation schemes. The wavelet transform results in a set of wavelet 
surfaces, one for each scale. We address the interception problem using two ap- 
proaches, either of them can be used independently, but they differ in performance as 
will be shown later. In the first approach, we visually inspect the 2-D wavelet surfaces 
to identify and classify the structure of the FH signal and obtain an estimate for the 
hop time interval. In the second approach, we apply a proposed processing scheme to 
estimate the hop start/stop times, the hop-scale pattern, and the hop frequency. The 
extraction of the hop start/stop times is addressed using an edge detection approach 
by applying a compass operator which is well known in the image processing area. 
The hop-scale pattern is obtained by applying an energy analysis. The energy anal- 
ysis assigns a scale index (called the proper scale) to each hop. The proper scale, for 
each hop, is that scale which has the greatest energy content. The sequence of proper 
scales, representing the hop sequence, is called the hop-scale pattern. The frequency 
of each hop can be extracted from the wavelet surface at the proper scale. 

The identification of the FH signal may be accomplished based on: 

(1) The hop-scale pattern; If two or more wavelet scales are applied, the hop-scale 
pattern of the FH signal will be different from the hop-scale patterns of other 
digital modulation signals. 

(2) The frequency diversity: If all frequencies reside in one scale, then, an FH sig- 
nal will have different frequency components as a function of the hop intervals. 

Recall that, the ICF of the FH signal has a cellular structure which comprises 
of a tiling of diamonds. Each hop results in a diamond with a width equal to the hop 
interval. Thus, the diamond boundaries point to the hop start/stop times. All wavelet 
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surfaces at all scales emphasize the signals and hence the diamonds which belong to 
that scale. Other signals are attenuated. In this chapter, we address aspects of the 
discrete implementation, the proposed processing scheme, the measured parameters, 
and performance measures. In section VLB, we address the discrete form of the 
IGF and its output expression to a white noise input. In section VI.C, we address 
the visual inspection and FH identification from the wavelet surfaces and obtain an 
estimate for hop time intervals. Then we compare the FH wavelet surfaces with 
wavelet surfaces from other digital modulation schemes. In section VI. D, we consider 
an energy analysis for identifying the scale of each frequency hop. We also investigate 
the performance of scale identification and evaluation measures. In section VI.E, 
we address the equalization of the spectral shape of the IGF and its impacts on the 
performance of scale identification. In section VI. F, we address the hop frequency 
extraction from the wavelet surfaces at the proper scale and from the original time 
signal. We also investigate the performance of frequency extraction and evaluation 
measures. Finally, we address the extraction of the hop start/stop times in section 
VI.G. 

B. PROCESSING SCHEME 

The interception problem usually assumes some prior knowledge about the sig- 
nal of interest. In our case, we assume the signal hopping bandwidth is approximately 
known and the data is properly heterodyned and sampled. 

1. Processing Sequence 

The input data parameters to our processing scheme are the parameters of the 
FH signal (i.e., sampling frequency and the range of possible frequencies) while the 
outputs are: 

(a) The wavelet transform of the IGF surface, i.e., the wavelet surfaces. 

(b) The scale of each hop. 

(c) The frequency of each hop. 
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(d) The hop interval start/stop times. 

The data (input) process generates FH signals according to some predeter- 
mined input parameters (i.e., number of hops, hop interval, hop frequencies, sampling 
frequency and SNR). The analysis process generates the ICF surface and computes 
its wavelet transform at the predetermined wavelet scales. The measurement process 
extracts the outputs (the scale of each hop, the frequency of each hop, and the hop 
interval start/stop times). 

Figure 14 shows the functional description of the processing scheme where the 
Input stands for the input process. The Hilbert Transform generates the analytic form 
of the signal. The Instantaneous Correlation Function generates the (discrete time) 
ICF surface. The Wavelet Transform computes the (discrete time) wavelet transform 
of the ICF surface (i.e., wavelet surfaces). The Scale Identifier identifies the hop scale. 
The Hop Frequency extracts the hop frequency from the wavelet surfaces. The Hop 
Timing extracts the hop start/stop times and the hop interval. 

2. Discrete-Time Implementation of the Instantaneous 
Correlation Function 

For discrete implementations, let R{n, u) define the ICF of the discrete-time 
signal x{n), 

R{n, u) = x(n + ^'^ x* ~ . (VI-1) 

where n is time and u is the time delay. Note that, the normalized sampling interval 
is 1 for the discrete implementation. In addition, the combined index n ± | should 
assume only integer values. Therefore, there will be two choices; 

• Zero-inserted ICF: If we let u assume only even integer values, then R{n,u), 

for odd u, must be set to zero. 

• Frequency-doubled ICF: If we let u = 2m then R{n,m) will be defined as 

R{n,m) = x{n + m)x*{n — m) . (VI-2) 
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Zero-inserted ICF. The effect of zero insertion to the ICF will result in doubling 
the number of wavelet coefficient over the resultant surfaces, implying more computa- 
tion and storage load. Each hop will reside within its scale according to its associated 
signal frequency. 

Frequency-doubled ICF. The frequency-doubled discrete ICF has the following 
characteristics in contrast to the zero-inserted ICF: 

• Due to doubling the frequency, the minimum sampling frequency of the mon- 
itored signal will be twice the Nyquist rate. 

• We obtain half the number of coeflBcients at the ICF surface and consequently 
at the wavelet surfaces. 

• Over the wavelet surfaces, each hop will be located at a lower scale index than 
its associated signal frequency (lower scale index indicates higher frequency 
according to the MATLAB designation, see VI.A.3). 

Since we usually pick a sampling frequency which is a multiple of the Nyquist 
rate, the choice of the frequency doubled ICF is advantageous because of a fewer 
number of coefficients, hence, saving in computation and storage. On the other hand, 
the loss of processing gain, due to pushing each hop to a lower scale index, can be 
overcome by chosing higher sampling frequency. For these reasons we adopted the 
definition of the frequency-doubled ICF in the processing scheme. 

3. MATLAB Discrete Wavelet Transform 

The MATLAB Wavelet Toolbox is used to implement the discrete wavelet 
analysis of the ICF surface. The wavelet transformation is carried out using the 
one-dimensional wavelet transform in the direction of the time delay u for each time 
element of the correlation function (as previously addressed in Chapter III). MATLAB 
uses a scale index designation and a corresponding multi-level decomposition tree 
as shown in Figure 15. Sj denotes the output of the wavelet transform at the 
wavelet scale (i.e., the wavelet coefficients or the details at the level). Aj denotes 
the approximation at the level. Hq is the scaling filter (low pass filter) and Hi is 
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the wavelet filter (high pass filter) [Ref. 59]. MATLAB uses the Daub-N designation 
where N indicates the wavelet order and the number of vanishing moments, while the 
actual filter length is 2 AT [Ref. 59). 

C. VISUAL IDENTIFICATION 

In this section, we address the visual inspection of the wavelet surfaces in 
order to identify the FH modulation scheme. We investigate different approaches to 
represent these surfaces and considered real and complex valued wavelet functions. 

1. Real and Complex- Valued Wavelets 

According to Daubechies [Ref. 42], there is no symmetric or anti-symmetric 
real compactly-supported orthonormal wavelet. Symmetry property can be achieved 
only for complex-valued wavelets. Symmetry of a wavelet implies that the FIR filter 
representation has a linear phase response [Ref. 36] . It is an important feature in some 
wavelet applications, such as numerical resolution of partial differential equations 
with boundary conditions [Ref. 60]. A complex- valued Daubechies wavelet of order 
3 is given in [Ref. 60] and has been used for comparison with the MATLAB real 
Daubechies wavelet of the same order. The Daub-3 complex wavelet has the following 
scaling coefficients: 

'.«(!) = M6) = 

M2) = A.(6) = 

ft„(3) = A„(4) = (VI.3) 

2. Surface Representation 

Structure of the FH signal can be identified and classified by visually inspecting 
the wavelet surfaces. We can also obtain an estimate for the hop start/stop times. The 
next automated step is achieved by applying a proposed processing scheme to estimate 
the hop start/stop times, the hop-scale pattern, and the hop frequency. Using the 
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analytic form of the FH signal, the ICF surface will be a complex-valued function. 
Therefore, for any real or complex wavelet function the resultant two-dimensional 
wavelet surfaces are complex-valued. Hence, we consider the real part, the imaginary 
part, the magnitude, or the phase (angle) of the complex-valued surface. 

An opinion test was carried out to assess the quality of identification from the 
wavelet surfaces using different representations. The Morlet wavelet and Daubechies 
wavelet of order 3 were both used in their real and complex form. Four different 
representations were investigated: real part, imaginary part, magnitude (absolute 
value), and phase. The objective of the opinion test is to identify the cellular structure 
over the wavelet surface and to identify the diamond boundaries as an indication of 
the hop start /stop times. 

Figures 16 to 19 present the real part of the WT obtained with the real-valued 
Daubechies wavelet of order 3, at an input SNR values of 10 and 3 dB. In Figures 
16 and 17 we note that we can not identify the FH structure over the surface of the 
ICF denoted by ”CF”. For wavelet surfaces, labeled ’’Sir”, ”S2r”,..., ”S5r”, we can 
identify diamond patterns at the hops number 1 , 2 ,..., 5 , respecively. Same findings 
can be observed in Figures 18 and 19 about wavelet surfaces. These results show that 
we can identify the FH structure from the wavelet surfaces, while it is not possible to 
do so from the ICF surface. 

D. ENERGY ANALYSIS AND SCALE 
IDENTIFICATION 

Energy analysis is performed over the wavelet surfaces at all scales considered. 
For the energy analysis we assume correct hop timing. The hop interval and hop 
start/stop times are estimated by visual inspection (and later on from the processing 
scheme). Thus, we can point to each hop over the wavelet surface and compute the 
energy contained. The energy contained at all scales are compared to identify the 
proper scale of each hop. 
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1. Energy and Energy per Sample 

Parseval’s theorem for the complete orthogonal filter bank (over L partitions) 
is applied to the discrete time wavelet analysis [Ref. 45]. It is given by: 

IW")lP= E (|C(L. 2 t)|' + EWj, 2 *: + l)N , (VI. 4 ) 

k&Z \ j=l ) 

where, in the sense of wavelet analysis, C{L,2k) are the scaling coefificients at the 
scale L, d{j, 2k + 1) are the wavelet coefficients at scale j, and k is the wavelet shift 
variable. 

Therefore, the quantity 

00 

m= E (vi-5) 

/:=—oo 

represents the portion of the signal’s energy over the scale. For narrowband signals 
which reside within one scale, the signal energy will belong to that scale. However, 
some portions of the signal energy will be leaking to other scales. To identify the 
proper scale one can consider the maximum value due to the following quantities: 



1) Wavelet coefficient. 

2) Total energy. 

3) Energy per sample. 



Energy per sample is defined as: 



AU) = 






(VI.6) 



where j is the scale index, E(j) is the total energy at the scale, and N(j) is the 
number of wavelet coefficients at this scale (i.e., the number of surface coefficients). 
According to the scale designation of the MATLAB Wavelet Toolbox, 



JVy) = + 1). 



(VI.7) 
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If the signal x{t) resides within the scale j or resides within the scale (j + 1), 
with the same total energy, we have 

A{j + 1) = 2A{j) . (VI.8) 

Thus, there will be a 3 dB gain in the energy per sample per each increment in the 
scale index. 

Table 1 summarizes the energy distribution obtained with wavelet analysis 
using Daubechies wavelet of order 2 (Daub- 2) and order 10 (Daub-10). There are 
three input vectors with 64 samples each, the first vector has a frequency of S/SF^, 
the second has 3/16Fs, and the third has 3 / 32 F 5 , where Fg is the sampling frequency. 
This means that the first input vector resides within the first scale, the second vector 
resides within the second scale, and the third vector resides within the third scale. 
We conclude the following: 

• The total energy of the input signals is distributed among the scales. The sum 
of the total energies over the scales is slightly less than the input signal total 
energy since we disregard contribution from the low pass section. 

• The proper scale (where the signal resides) has the greatest share of the total 
energy. This share is greater (percentage wise) if a longer wavelet (e.g., Daub- 
10) is used than for a shorter wavelet (e.g., Daub-2). 

• Energy per sample at the proper scale is larger if the signal resides at a higher 
scale. The gain factor in energy per sample (if the signal resides within scale 2 
rather than scale 1) is approximately 1.41 and 1.64 for Daub-2 and Daub-10, 
respectively. We note that ideally, the gain factor in energy per sample should 
be 2 per scale index. 



72 



Table 1: Energy Distribution of Daub-2 and Daub-10 Wavelets 



Measure 


Wavelet 


Scale 


/ = (3/8)F, 


Signal with 
/ = (3/16)F, 


/ = (3/32)F, 






1 


1.3365 


0.6574 


0.4349 




Daub-2 


2 


0.5753 


1.6998 


1.0356 


Max 




3 


0.3743 


0.3594 


1.8480 


Coefficient 




1 


1.2524 


0.3818 


0.3284 




Daub-10 


2 


0.3968 


2.0331 


0.5504 






3 


0.1854 


0.3859 


2.8296 






1 


29.7687 


7.2565 


0.8486 




Daub-2 


2 


1.1313 


23.3282 


8.2775 


Total 




3 


1.0495 


0.9492 


21.2601 


Energy 




1 


31.4508 


1.6228 


0.4040 




Daub-10 


2 


0.3980 


29.9017 


1.9671 






3 


0.0875 


0.3900 


29.0447 






1 


0.9201 


0.2199 


0.0257 


Sample 


Daub-2 


2 


0.0628 


1.2960 


0.4599 


Energy 




3 


0.1049 


0.0949 


2.1260 


(average 




1 


0.7671 


0.0396 


0.0099 


per sample) 


Daub-10 


2 


0.0133 


0.9967 


0.0656 






3 


0.0036 


0.0163 


1.2102 
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2. Scale Identification in an Ideal Noise Environment 
In this particular simulation, the correlation surface is approximated by si- 
nusoidal signals embedded in ideal white Gaussian noise. The maximum coefficient, 
total energy, and energy per sample are tested. The signal-to-noise ratio is expressed 
as: 

_ input signal power 
input noise variance 

Three different input vectors at SFg/S, 3Fs/16, and 3Fg/32, each Avith 2048 samples, 
are analyzed using Daub-2 and Daub-10 wavelets. SNR values are 10, 5, 2.5, 0, —2.5, 
—5, —10, and —15 dB. Each SNR level was used in 20 trials. Results are shown in 
Table 2 where the tabulated SNR is the smallest SNR at which each measure correctly 
designates the proper scale 100% of the time. 



Table 2: Noisy Observation 



Measure 

Wavelet 

Type 


Greatest 

Coefficient 


Total 

Energy 


Sample 

Energy 


Smallest SNR [dB] 


Daub-2 


2.5 


0 


-10 


Daub- 10 


0 


-2.5 


-10 



These results show that the energy per sample achieves the lowest SNR value com- 
pared to the other two measures. Consequently, this measure is considered the most 
reliable one for proper scale identification. 

3. Hop-Scale Pattern 

We recall that scale identification means assigning a scale index to each hop 
of the observed signal. The hop is assigned the scale whose energy per sample is the 
greatest among the values of the other scales; hence, it is called the proper scale. 
Therefore, a sequence of hops will have a sequence of proper scales, also called the 
hop-scale pattern. 
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Properly selecting a heterodyne and a sampling frequency will result in proper 
choice of the wavelet scales. The number of scales is defined by the number of octaves 
in the heterodyned signal. The designation of the scales (scale indices) is controlled 
by the selected sampling frequency, as each scale index is defined in terms of fractions 
of the sampling frequency as shown in Figure 15. As an example, assume we have 
an FH signal with hopping bandwidth of 30-90 MHz. If the signal is heterodyned to 
10-70 MHz and sampled at 280 MHz, then the proper wavelet scales are S2, S3, and 
S4. If the FH signal is heterodyned to 2-62 MHz and sampled at 280 MHz, then the 
proper wavelet scales are S2, S3, S4, S5, and S6. 

To distinguish the FH signal from other modulation types, the number of 
scales used must be at least two. The hop-scale pattern will then be comprised of 
two symbols (i.e., two scale numbers). For other modulations the hop-scale pattern 
will be comprised of one symbol (i.e., one scale number). For the MFSK scheme with 
a carrier frequency in the above mentioned band, the bandwidth of the modulated 
signal is typically 25 KHz (i.e., the channel separation for the wireless sets). Thus, 
the MFSK signal bandwidth can not span an octave except when we heterodyne to 
a very low frequency (i.e., less than 25 KHz). Hop-scale patterns are also useful for 
hop frequency extraction from wavelet surfaces, since they point to the proper scale 
where most of the signal energy is concentrated. Therefore, frequency extraction, if 
done on that wavelet surface, will provide correct values. 



4. Success Rate 

Performance of scale identification is evaluated as the success rate, Pid- We 
generate known hop-scale patterns and obtain the percentage of the correctly identi- 
fied hop-scales. Therefore, Pid is defined as: 



Pid = 



number of correct hop-scales 
total number of hops 



(VI.9) 



The probability of error, given by Pe = 1 — Pj^, is the error in identifying the correct 
scale. The quality of scale identification depends on the height of the greatest energy 
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per sample relative to the other sample energies from other scales. Therefore, a 
suggested measure, dh, of the quality of scale identification is defined as: 



^ max({^(j)}) - mean({^(i)}) 

\/var({^(i)}) 

where max(-) is the maximum value, mean (•) is the mean value, and var(-) is the 
variance. The term d^, measures the distance between the height of the energy per 
sample at the proper scale to its average values over all scales, and is expressed in 
terms of the standard deviation of the energy per sample. 



(VI.IO) 



E. EQUALIZATION OF THE SPECTRAL SHAPE OF 
THE INSTANTANEOUS CORRELATION FUNCTION 

Simulations show that the IGF of white noise (data set) has a triangular type 
spectrum in the direction of the time delay. Figure 20 shows theoretical and experi- 
mental frequency spectra taken in direction of the time delay (i.e., transforming the 
time delay to frequency). Thus, the energy per sample for all wavelet scales need to 
be compensated since the spectrum is colored. 

Figure 20 shows that the slope of the IGF spectrum is 1, thus, its height, q{j), 
at the middle of the j*** scale is: 



qU) = 1 + 



3(iV - 1) 

23+1 ’ 



where N is the number of IGF data points in the direction of the time delay. Note 
that the IGF of white noise is created by multiplying samples of the white noise in 
the time domain. Therefore, the resultant FT of the IGF is the convolution of two 
FT of the white noise, consequently, we obtain the triangular spectral (FT) shape. 
As a result, we will apply an equalizer, whose values are the reciprocal of q{j), to 
compensate the magnitude distortion at all scales. 

The performance of the scale identification is evaluated and a new success rate 
Pid is computed. Results are presented Ghapter VII. 
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F. FREQUENCY ESTIMATION 

In Chapter III, we related the WVD of the signal x{t) to the wavelet trans- 
form of its ICF. The relation is stated in Equation III.33 and it is repeated here for 
convenience: 

WUt-,sJ) = r 

J—oo 

where is the wavelet transform of the ICF of x{t), Vx{t,f) is the WVD 

of x{t). In Chapter 2, we discussed the properties of WVD and concluded that the 
WVD provides the means to estimate the signal frequencies. The FT of the wavelet 
surface gives a bandpass filtered version of the WVD of the FH signal. Thus, the hop 
frequency can be extracted from a Fourier transform of the wavelet surface, in the 
delay direction, over the main diagonal of the proper diamond. 

1. Frequency Resolution 

The Fourier transform of the wavelet coefficients can be used for spectral 
estimation at any scale. Consequently, the frequency resolution will be dependent 
on the Fourier transformation. The frequency resolution (i.e., the minimum spacing 
between two resolved narrow band components) of the DFT is approximately equal 
to A/ = Fg/N. Thus, at any given scale k, the number of data points N{k) will be 
related to the number of input data points N by: 

m) = -wi. 

The sampling frequency at which the data points are taken is scale dependent, i.e., 
Fg(k) = Fg/2^, where Fg is the input sampling frequency. 

At the scale, both the number of data points (wavelet coefficients) and the 
sampling frequency have been reduced by the same factor. Consequently, the FT of 
the N{k) data points has a frequency resolution given by: 

This means the frequency resolution of the Fourier transform is constant, independent 
of the wavelet scale. 
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2. Success Rate 

Performance of the frequency estimation procedure is evaluated in terms of the 
success rate, Pf. The hop frequency is considered correctly estimated if the spectral 
peak is at the true spectral bin. Pf is defined as: 

number of correct hop frequencies 



Pf = 



(VI. 12) 



“'r 

total number of hops 

The probability of error in estimating the correct hop frequency is denoted by Pg, 
where Pg = 1 — P/- The hop frequency is estimated as the bin number corresponding 
to the peak of the FT over a specified region of the wavelet surface. Therefore, the 
quality of the frequency estimation depends on the spectral height of the peak relative 
to the average background. A measure of the quality of the frequency estimation is 
given by: 



df = 



max(|X(/)|) - mean(|Ar(/)|) 



(VI. 13) 



where X{f) is the FT over the wavelet surface, max (•) is the peak value, mean (•) is 
the mean, and var (•) is the variance. The variable df measures the distance between 
the peak value of X (/) to the average background in units of the standard deviation 
of X if). 



3. Alternatives for Frequency Estimation 

Given the correct hop start/stop times, the hop frequency may also be esti- 
mated directly from the time signal or from the IGF surface. To extract the frequency 
from the original time signal we use the FFT over the hop length. Recall that the 
FFT is a matched filter for sinusoidal signals in white Gaussian noise. Thus, a bet- 
ter performance is expected in comparison to the nonlinear processing of the signal 
through the IGF computation and the wavelet transformation. The performance of 
frequency estimation from the signal or from its IGF is evaluated using the already 
defined measures P/ and df. We note that Pfs and Pfc denote the success rate using 
the original time signal or using the IGF surface, respectively. Also, dfs and dfc de- 
note the measure of the frequency estimation quality at the original signal or at the 
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ICF surface, respectively. The performance of Pfs and P/c are evaluated and results 
will be given in Chapter VII. 

G. EXTRACTION OF HOP TIMES 

Hop timing extraction is obtained by estimating the hop start and stop points. 
We recall that ICF surface and wavelet surfaces have both a cellular structure consist- 
ing of diamonds, where each diamond is associated with a specific hop. The diamond 
edges define the start /stop point of the hops. The diamond width corresponds to 
the hop interval Th. The sides (edges) of the .diamonds of hops are mutually par- 
allel and spaced by the hop interval Th- Therefore, we obtain the distance between 
two sequential intersections of the diamonds and the time axis, to determine the hop 
start (or stop) point. This tends to transform the transition point detection into an 
edge detection problem. There are many approaches to tackle the problem of hop 
timing extraction; in the following section we address one technique based on an edge 
detection operator. 

Compass Operator 

Edge detection is a fundamental problem in image analysis since edges help 
to identify objects. Edges are characterized by an abrupt change in gray level, there- 
fore, edges can be detected using the derivative (gradient) operators which will be 
maximum in the direction of the edges. There are two types of edge operators, the 
gradient operators and the compass gradient operators (also called the compass oper- 
ators) [Ref. 62]. The gradient operator measures the gradient of the two-dimensional 
image in two orthogonal directions. It is usually applied to detect edges with unknowm 
directions. The compass operator measures the gradient of the two-dimensional im- 
age in a specific direction. It is used to detect edges with pre-determined directions. 
Since the wavelet surface has edges at angles of ^ and ^ radians, a specific com- 
pass operator can be used to detect the edge in these directions. There is a variety 
of compass operators. They differ in form, depending on the direction of the edge 
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to be found. An example of the compass operator is the NE compass ( NE stands 
for North-East). It has a 3 x 3 matrix NE = [0, 1, 1; — 1,0, 1; — 1, — 1, 0]. Another 
example is the NW compass (NW stands for North-West) with the 3x3 matrix 
NW = [1, 1, 0; 1, 0, —1, 0, —1, —1]. The compass operator is applied to the upper half 
of the wavelet surfaces to potentially detect the diamond edges. An array of the NE 
compass operator is shown in Figure 21. The compass array is used to scan the sur- 
face from left to right. All of the contributions are summed according to the weights 
of the 3x3 matrix. The maximum value is extracted to define the point where the 
compass array matches an edge. To make our data applicable to compass operations 
we need to take care of the negative portions of the surfaces, which is done by adding 
in a positive number equals in magnitude to the smallest (i.e., the most negative) 
surface value. Results are presented in Chapter VII. 
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Figure 14. Functional description of the processing scheme. 
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Figure 15. MATLAB wavelet scale designation and computation tree. 
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Figure 16. Wavelet surfaces of FH signals using Daub-3 at 10 dB (scale 1, 2, and 3). 
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Figure 17. Wavelet surfaces of FH signals using Daub-3 at 10 dB (scale 4 and 5). 
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Figure 18. Wavelet surfaces of FH signals using Daub-3 at 3 dB (scale 1,2, and 3). 
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Figure 19. Wavelet surfaces of FH signals using Daub-3 at 3 dB (scale 4 and 5). 
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Figure 20. Theoretical and simulated spectrum of the IGF of white noise. 
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Figure 21. The compass line over the upper half of the wavelet surface. 
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VII. 



SIMULATIONS AND RESULTS 



A. INTRODUCTION 

In this chapter, we will provide experimental results of the processing scheme 
introduced in Chapter VI. Section VII.B shows the results obtained by visual inspec- 
tion of wavelet surfaces as they pertain to the identification of FH signals. Wavelet 
surfaces from other digital modulation schemes are compared with those resulting 
from the FH signals. Section VII.C provides results of hop-scale identification. Sec- 
tion VII. D provides results of hop frequency extraction and compares results obtained 
using frequency estimation from wavelet surfaces and from the original time signal. 
Section VII. E provides the results of hop start/stop times estimation obtained by 
using the compass operator over the wavelet surface. 

B. VISUAL INSPECTION 

We carried out an opinion test to identify the FH scheme by examining the 
wavelet surfaces. Ten participants were involved, each one was asked to identify the 
diamond structure of the FH signal over the wavelet surfaces at all pertinent scales and 
for all hops. Two types of wavelet were used; the Morlet wavelet and the Daubechies 
wavelet of order 3. Both wavelets were used in their real and complex form. Four 
SNR values were used; 10, 6, 3 and 0 dB. Four different surface representations 
were used; the real part, imaginary part, magnitude, and phase. To avoid biasing 
and preconception all participants started to identify the surfaces resulting from the 
lowest SNR value and then the higher SNR values consequently. Detailed scoring 
tables and the scoring code are given in Appendix B. Only scoring tables belonging 
to 10 and 3 dB are given in this section because of their significance to the final 
conclusion of the opinion test. 10 dB is the highest tested SNR value while 3 dB 
is the minimum SNR which provided the minimum acceptable identification score of 
2 (according to the scoring code). In Appendix B, each scoring table is dedicated 
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for a specific wavelet type at a specific SNR value and contains the scores given for 
different representations at different wavelet scales. ’’Real” stands for representing 
the wavelet surface by the contour plot of its real part, while ”Imag.” , ”Abs.” and 
’’Angle” stand for representing the surface using the contour plots of imaginary part, 
magnitude and phase, respectively. The fractions in numbers listed in the tables came 
from averaging the ten scores given by the ten participants. 

Table 3 shows the summary of scores obtained from different wavelet types at 
SNR values of 10 and 3 dB when representing the wavelet surfaces by their real part. 
Eventually, the real part and the imaginary part outperformed other representations. 
The ratings in Table 3 ranges from 0.2 to 1, where 1 represents perfect identification 
of the hop diamonds at their proper positions and 0.2 represents just distinction of 
the hop diamonds from the background noise. Test signals were designed to reside 
within the five wavelet scales given in the table. Comparing different scoring tables 
in Appendix B, it becomes apparent that the real part or imaginary part provided 
the best surface representation for the purpose of visual inspection. The magnitude 
or phase provided a very poor representation. We noted that real valued wavelets did 
a slightly better job than the complex valued wavelets (at least for the wavelet types 
used in the simulations). 

Thus, this visual opinion test indicates that: 

(1) The FH signal can be identified over the wavelet surfaces by its cellular struc- 
ture which is dominant at the proper scale. That is, each scale will emphasize 
the hops which belong to the scale and attenuate other hops. The (diamond) 
cellular structure provides an estimate for the hop start/stop times. 

(2) The FH signal can be well identified at 3 dB SNR and above. 

(3) The real part or imaginary part provides the best surface representation for 
the purpose of visual inspection. 

(4) The real form of the wavelet function provides better surface representation 
than the complex form for the wavelets used in simulations. However, other 
types of wavelets may perform differently. 
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(5) The Morlet wavelet has better performance than the Daubechies wavelet es- 
pecially at higher scales. This is because our Morlet wavelet has a narrower 
bandwidth than that obtained with Daubechies wavelet, and the Daubechies 
wavelet results in a decimated transform which reduces the number of coeffi- 
cients at higher scales. 

(6) Other modulation schemes such as ASK, PSK, and MFSK will have wavelet 
surfaces residing only at one scale. Plots of these wavelet surfaces are given in 
Appendix C. 
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Table 3: Summary of identification scores for real part representation. 



Wavelet 


SNR 


Scales 


Type 


[dB] 


SI 


S2 


S3 


S4 


S5 


Real 


10 


1.0 


1.0 


1.0 


0.95 


0.95 


Morlet 


3 


1.0 


1.0 


1.0 


0.95 


0.9 


Complex 


10 


1.0 


1.0 


1.0 


0.9 


0.75 


Morlet 


3 


0.7 


0.8 


0.9 


0.5 


0.4 


Real 


10 


1.0 


1.0 


0.8 


0.2 


0.2 


Daub-3 


3 


0.8 


0.7 


0.5 


0.2 


0.2 


Complex 


10 


1.0 


1.0 


1.0 


0.3 


0.2 


Daub-3 


3 


0.6 


0.5 


0.7 


0.4 


0.2 



• Scores vary between 0.2 and 1, where 

• 1: perfect identification of hop diamonds. 

• 0.2: just distinction of hops from background noise. 
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C. SCALE IDENTIFICATION 

The proposed processing scheme is used to extract hop start/stop times, the 
hop-scale pattern, and the hop frequency. First, we investigate scale identification 
and hop frequency extraction assuming correct hop timing is available. Extraction 
of hop start/stop times will be considered later. For scale identification performance 
we carried out simulations using the following data: 

(1) Signal pattern length: 5 hops. 

(2) Wavelet scales: SI, S2, S3, and S4. 

(3) Wavelet types: Daub-2, 4, and 8. 

(4) SNR: from 10 to —10 dB. 

(5) Number of experiments: 250 per wavelet per SNR value. 

The hop frequencies of FH signals are selected to generate the hops according 
to a known scale pattern taking in consideration the doubling effect of frequencies by 
computing the ICF surface. The test pattern is selected of four scales; SI, S2, S3 and 
S4. The wavelet surfaces are generated from the ICF surfaces at those relevant scales. 
The total energy of each hop at each scales are computed and the energy per samples 
are obtained by dividing the total energy by the number of wavelet coefficients at 
each scale. Over each hop, the scale with the greatest energy per sample is estimated 
as the proper scale. The resultant estimated hop pattern is compared to the known 
hop pattern and the probability of correct identification is computed. We compute 
the energy contained in each hop over the diamond corresponding to that hop. We 
disregard energy contained outside the diamond since it is contributed by different 
hops. To avoid bias from the colored noise of the ICF surfaces an equalization must 
be performed to the resultant energy per samples at all pertinent scales before using 
them in deciding upon the proper hop scale. Although, the equalization process is 
based on the theoretical spectral shape given in Figure 20, the resultant equalized 
energy per samples, for ideal input white Gaussian noise, do not have the same value 
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at all scales. Table 4 presents averages and variances obtained for the equalized 
energy per samples for Daubechies of order 2, 4, and 8. We note that the average 
values for Daub-2 and Daub-8 are monotonically increasing with the scale index. For 
Daub-4 they are monotonically increasing for the scales S1,S2 S3 and S5 while they 
have a dip at scale S4. These different trends suggest that different wavelets respond 
differently to the theoretical equalization of the spectral shape of the ICF. Therefore 
the correction weights of the energy per sample depend on the wavelet used. In these 
simulations, results are given for the proper scale identification after equalization and 
correction. 

Figures 22 to 24 show the success rate, Pid obtained for Daubechies wavelets 
of order 2, 4, and 8 at scales SI, S2, S3, and S4. For the performance of scale 
identification we consider the minimum SNR at which Pid is 1 as the figure of merit. 
Over all tested scales the success rate, Pid assumes the value of 1 at different minimum 
input SNR values. This is a function of the order of the wavelet and the scale. Figure 
22 shows that the performance of Pu obtained from Daub-2 has achieved the value 
of 1 at SNR equals —1 dB at most of the scales, hence , —1 dB is considered the 
minimum SNR value for Daub-2. The minimum SNR value for Daub-4 and Daub-8 
is —2 dB at most of the scales as shown in Figures 23 and 24. We note also that 
for Daub-8, Pid assumes the value of 0.9 for most of the scales at an SNR of —3 dB. 
Results show that: 

(1) For long wavelets (Daub-8) the probability of correct scale identification is 1.0 
at SNR equals to —2 dB, while it is 0.9 at SNR equals —3 dB. 

(2) Longer wavelets perform better than shorter ones. 

(3) The performance at higher scales is mostly better than the performance at 
lower scales in terms of minimum SNR at which Pid is 1. The exception of 
that case is the performance at the scale Si which can be justified by the 
imperfect equalization of the ICF spectral shape. 
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Table 4: Mean and standard deviation of energy per sample. 



Mean 


SI 


S2 


S3 


S4 


S5 


Daub-2 


1.5688 


1.7557 


1.9479 


2.3165 


2.6086 


Daub-4 


1.6303 


1.8432 


2.0218 


2.5585 


2.1149 


Daub-8 


1.7173 


1.9979 


2.2725 


2.4776 


2.7525 



Standard 

Deviation 


SI 


S2 


S3 


S4 


S5 


Daub-2 


0.3419 


0.3763 


0.6503 


0.9542 


1.3176 


Daub-4 


0.3507 


0.4103 


0.6891 


1.1328 


1.2776 


Daub-8 


0.3665 


0.4520 


0.7874 


1.1063 


1.5845 



Theoretically, wavelet surfaces at higher scales have higher SNR values than 
those at lower scales due to the higher processing gain at higher scales. Consequently, 
the decay in the Pid performance, at low input SNR, is slower at higher scales than 
the decay at lower scales. But, we noted that the decay in Fid is not consistent 
at different scales, especially at scale SI, which shows, mostly, the slowest decay 
compared to other scales. This anomaly may come from many interacting effects due 
to the imperfect equalization and the non-ideal wavelet filters at different scales. 

D. FREQUENCY EXTRACTION 

We carried out simulations to evaluate the performance of frequency estimation 
using the data from Section VII. C. 
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The hop frequencies are estimated by taking the FT of the vector of wavelet 
coefficients located at the center of the diamond in the direction of the time delay. 
The bin corresponding to the peak value represents the estimated hop frequency. The 
estimated hop frequency is then compared to the true hop frequency and the proba- 
bility of correct frequency extraction (the success rate) is computed. The estimated 
frequency is considered correct if the estimation error is less, in percentage of the true 
frequency, than where N is the length of the vector of the wavelet coefficients. 

Figures 25 to 27 plot the success rate P/ obtained for different wavelets at 
different scales. Figures 28 to 30 display the corresponding values of the distance df. 

We consider the minimum SNR, at which P/ is 1, as the figure of merit for 
the performance of frequency extraction. Over all tested scales the success rate of 
frequency estimation, P/, assumes the value of 1 at different minimum input SNR 
values. 

Figure 25 shows that the P/ value for Daub-2 is 1 at SNR equals 0 dB at most 
of the scales, hence, 0 dB is considered the minimum SNR value for Daub-2. The 
minimum SNR value for Daub-4 and Daub-8 is also 0 dB at most of the scales, as 
shown in Figures 26 and 27. 

For Daub-2, Daub-4 and Daub-8 the minimum input SNR is 0 dB at scales 
Si, S2, and S3. At a success rate P/ of 0.9, the minimum SNR equals —2 dB over 
most scales. The distance df is about 4 times the standard deviation of the wavelet 
surface data for all tested scales as shown in Figures 28 to 30. Values of the success 
rate P/ and the distance df shows that the hop frequency can be reliably extracted 
from the wavelet surfaces at an SNR equals 0 dB or better using only one FT at the 
center of the hop diamond area in the direction of the time delay. 

Alternatives for frequency estimation 

Hopping frequencies may also be estimated from the original signal or from the ICF. 
Figure 31 plots the performance Pfs and Pfc obtained from the time signal and from 
the ICF, respectively, at different SNR values, assuming exact estimates of hopping 
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start/stop times are available. Figure 32 plots the distance dfs and dfc obtained 
from the time signal and from the ICF, respectively, at different SNR values. The 
frequency estimation success rate from the time signal, Pjs, is 1 at SNR value of 
—5 dB. The frequency estimation success rate from the ICF, Pfc, is 1 at SNR value 
of 0 dB. The distances d/s and dfc, obtained from the original signal and from the 
ICF are 6 and 7, respectively. This shows that, assuming exact estimates of hopping 
start/stop times are available, hop frequencies may be estimated by processing the 
original signal at lower SNR values than can be achieved using the wavelet surfaces. 
The benefit obtained by analyzing the ICF surface by wavelet analysis is significant 
in case of unknown hop start/stop times. Also, less computations are needed, in case 
of estimating the frequency from wavelet surfaces, due to applying one FT per hop, 
using few coefficients. 

E. EXTRACTION OF HOP TIMES 

This section presents hopping time extraction results obtained using the com- 
pass operator discussed earlier in section VI.G. The wavelet surface is represented 
by its upper half plane resulting in a triangular shape instead of the familiar hop 
diamond shape. The line compass operator is used over the surface from left to right 
and the peak value of the resultant contribution is located over the time axis. The 
location of the peak value points to the estimated hop starting time. The difference 
between the true and the estimated starting time is considered the estimation error 
and is evaluated in terms of points over the time axis. Each SNR value is used in 20 
trials. Figures 33-35 show the results of the estimation error (Er) at different SNR 
values for different Daubechies wavelet lengthes. The mean-square error (MSE) and 
standard deviation (SD) of the estimation error are also given. Results show that 
the hop start/stop times can be extracted within accuracy of ±7 points (out of the 
128-point hop interval) at SNR values of 6 dB or better. 



97 






Figure 22. Pid for Daub-2. 
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Figure 23. Pid for Daub-4. 
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Figure 24. Pid for Daub-8. 
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Figure 25. Pf for Daub-2. 
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Figure 26. Pf for Daub-4. 






Figure 27. Pf for Daub-8. 
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Figure 28. The spectral distance df for Daub-2. 
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Figure 29. The spectral distance df for Daub-4. , 
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Figure 30. The spectral distance df for Daub-8. 
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Figure 31. P/s using the signal and P/c using the ICF. 
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Figure 32. dfg using the signal and dfc using the ICF- 
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Figure 33. Hop timing estimation error for Daub-2.^ 
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Figure 34. Hop timing estimation error for Daub-4, 
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Figure 35. Hop timing estimation error for Daub-8^ 
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VIII. CONCLUSIONS AND 
RECOMMENDATIONS FOR FUTURE WORK 

A. CONCLUSIONS 

Wavelet analysis of correlation functions is a new area which can be used 
to intercept communication signals embedded in noise. This work aimed at apply- 
ing wavelet analysis to the instantaneous correlation function to identify frequency 
hopped signals. 

The pure FH is the most popular scheme for frequency hopped signals. The 
instantaneous correlation function (ICF) of the complex-valued pure FH signal was 
shown to have a cellular (diamond) structure, where each hop contributes to one main 
diamond. Inside this diamond, the ICF has a single complex exponential component 
representing the hop frequency in the delay direction. We showed that the diamond 
intersections with the time axis are the hop start/stop times while the width of the 
diamond is the hop interval. 

The wavelet transform of the ICF surface results into a number of surfaces 
each at a specific wavelet scale (called the wavelet surface). The wavelet surface at 
any scale emphasizes the hops which reside in it and attenuates other hops. 

We addressed the interception problem using two possible approaches. In the 
first approach, we visually inspected the wavelet surfaces to identify and classify the 
structure of the FH signal and to obtain a rough estimate for the hop time interval. 
In the second approach, we applied the processing scheme which can also be used to 
automate the interception task. The proposed processing scheme estimates the hop 
start/stop times, the hop-scale pattern, and the hop frequency. The extraction of the 
hop start/stop times was addressed using an edge detection approach, by applying 
a compass operator, which is a well known technique in the image processing area. 
The hop-scale pattern is obtained by applying an energy analysis. 

The frequency of each hop can be extracted from the wavelet surface at the 
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proper scale, or as an alternative, from the original time data using the time param- 
eters. 

Visual inspection of the wavelet surfaces showed that the FH signal can be 
identified from the wavelet surfaces for an input SNR of 3 dB and above. Other 
modulation schemes such as ASK, PSK, and MFSK will have significant wavelet 
surfaces at one scale only since their frequency bandwidth does not span more than 
one wavelet scale. Hop timing extraction showed that the hop start/stop times can 
be extracted within accuracy of ±5% for short wavelets (Daub-2) at SNR of 6 dB or 
better. 

For the hop-scale identification, we adopted the energy per sample as a measure 
for the proper scale, assuming exact estimates of hop start/stop times are available. 
Results showed that the probability of correct scale identification is 1.0 at an input 
SNR of —2 dB for long wavelets (Daub-8), and it is 0.90 at an input SNR of —3 
dB. The performance of longer wavelets is better than that of shorter ones since 
longer wavelets have better spectral energy concentration than shorter ones. For 
hop frequency extraction, the success rate of frequency estimation from the wavelet 
surfaces showed that the probability of correct frequency extraction is 1.0 at an input 
SNR of 0 dB and above. 

Results showed that the minimum SNR required for automating the intercep- 
tion task is 6 dB. However, a better performance of frequency estimation is obtained 
at lower SNR value (—5 dB) by processing the original time signal, using the hop 
time information. 

B. RECOMMENDATIONS FOR FUTURE WORK 

For future extension of this work we recommend the following: 

(1) Addressing the automatic recognition of the cellular structure of the FH signal 
over the wavelet surfaces. There are two topics in the image processing area 
which may be helpful in this task; the automatic recognition of regions of 
interest, and automatic target recognition [Ref. 63]. 
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(2) Improving the performance of the hop-scale identification at lower SNR levels. 
In addition, it is suggested to readdress the equalization of the spectrum of 
the ICF surfaces. 

(3) Investigating other wavelet types, and the use of other definitions for the in- 
stantaneous correlation function. 

(4) Combining information from different wavelet surfaces to improve parameter 
estimation. 
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APPENDIX A. THE INSTANTANEOUS 
CORRELATION FUNCTION (ANALYTIC 

APPROACH) 



Let the complex FH signal be defined as: 

N 

E 

n=l 



s(t) = exp J + On) 






(A.l) 



where is the hopping frequency (in radians) and On is the carrier phase chip of the 
hop. N is the total number of observed hops. The shorthand notation t € {n} 
stands for the condition: 



to + {n- l)Th <t <to + nTh , 



(A.2) 



i.e., it confines the value of the time t corresponding to the hop duration. Th 
denotes the hop interval and to denotes an initial misalignment. The resulting in- 
stantaneous correlation function(ICF) is defined in terms of the signal values at time 
(t + t/2) and (t — r/2). 

If “t” is restricted over the hop, the sum term (t ±r/2) may stay over the 
hop or may extend to (n ± 1) hop. Consider a positive “r”, i.e., 0 < r < T/i, 
and let t be restricted to the first half of the hop and let to = 0) then, t will be 
restricted to 

(n-l)T;,<t < (A.3) 

and this is denoted as: 

t G , (A.4) 



and implies that if t is restricted to the 1®* half of the hop s{t -1- r/2) will remain 
within the hop since 0 < r < T^. Therefore, 



(t + |)=exp(i|;|c.„(( + I) 



t+(T/2)6{n), t€{(l/2)n} , 



(A.5) 
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If [t — Tj2) stays within the 'n}^ hop or the (n — 1)®’^ hop, then, 



+w„_x (t-|) 



t-(T/2)€{n}, te{(l/2)n} 

I) • (A.6) 

t-(T/2)€{n-l}, t€{(l/2)n}jy 

The condition t € {l/2n} is common for Equations A.5 and A.6. Therefore we take 
it from the equation body and place it aside as follows: 

A 



R‘((.r) 



expj ( (t + ^) 

\n=l I ^ 2/ t+(r/2)€{n} \ 2/ lt-(T 



/2)€{n> 



0 - i) 



t 6 



(H- 



(A.7) 

i— (r/2)€{n-l} ) J L ^ J 

The resultant terms of this summation depend on the validity of the indica- 
tor functions, i.e., terms to be combined together must have common (intersection) 
regions of their indicator functions. Thus, 



= expj (t + l) -w„(t-0 

(< + !)- 5) 



t+(r/2)6{n}, t-(T/2)e{n} 



t € 



(H- 



t+(T/2)€{7i}) t— ( t/2)€{ti— 1} 



1)^ 



(A.8) 



Let us examine the regions satisfying the following conditions: 

T T 

Condition 1 : t + — £ {n}, t — — £ {n}, 

Condition 2 : t+^€ {n}, ^ ^ € {n — 1}, 



for t € 






,T > 0. 



(A.9) 



The region where Condition 1 is valid can be defined from Equations A. 2 and 
A. 10 for {to = 0) as: 

{n-l)Th<t+'^<Th and (n - 1)T,, < t - ^ < nT^ , (A.IO) 
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where 



(n - \)Th <t< and Th> r >0. (A. 11) 

Equation A. 10 and A. 11 can be plotted according to the following set of inequalities 
(take, e.g., n = 1): 

^ + 2 ^ '^h, i + 2 — 

t > 0 , t < Q) n, 

0< r<Th. (A.12) 



An illustration of the corresponding set of equalities is shown in Figure 36. 
Therefore, condition 1 is valid over the triangle ABC. Similarly, condition 2 is 
defined by the following set of inequalities: 



t + ^<Th, t + ^>0, (A.13) 



and 0 < t{Th/2), 0 <r <Th. Condition 2 is valid over the triangle ADC, therefore 
t) = expj ujn (^ + 0 - (t - 0 



Condition 1 



\ \ Condition 2/ 



or 



I ^ 

R\t,r) = expj ( ^o;„r 

n=l 



+ 



Condition 1 



{Uji — U}n-i)t + (oJn + ^n-l) 



2 . Condition 2 J 

(A.14) 



where 0 < r, t € {(l/2)n}. 

Note for all hops, condition 1 and condition 2 valid areas will be as shown in Fig. 37 
where the areas of valid condition 1 are denoted by 1 and the areas of valid condition 
2 are denoted by 2. 
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For the case of negativ r, while t is still in the 1®* half of the hop, we make 
use of the conjugate symmetry property of the ICF, this will conclude Condition 3 
and Condition 4 as shown in Figure 38. When t is within the second half of the 
hop, similar procedure will lead to the other 4 conditions ( i.e.. Condition 5, 6, 7, 8). 
These areas are shown in Figure 39. 

Let us now summarize our observations from the derived formulai of the ICF 
by introducing the doubly indexed function Rm,n(t,T)- will represent the 

ICF of the FH signal over the t — r plane inside and outside the diamond of each hop. 
Rm,nit,T) will give the ICF expression once we assign its indices (m,n). The indices 
are shown in Figure 9. t) is given by 

^,n(i,r) = expi|(w,„ (WTO + Wn)!! , (A.15) 

where m and n are the indices of the two hops. 
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X 




Figure 36. Geometric definition of Condition 1 and Condition 2. 
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Figure 37. Areas of Condition 1 and Condition 2. 
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Figure 38. Areas of Condition 3 and Condition 4. 
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Figure 39. Areas of Condition 1 , 2 ,..., 8 . 
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APPENDIX B. SCORES OF OPINION TEST 



This appendix includes the tables of scores of the opinion test carried out for 
visual inspection and identification of the FH signal. Only one table is shown for each 
wavelet type. Each table is the average of 10 corresponding tables originated from 10 
experiments. This justifies the fractions in the tabulated scores. 

Table B.l: Wavelet: Complex Daub-3, SNR:10 dB 



s 


Real 


Imag. 


Angle 


Abs. 


1 


5 


5 


4.25 


1.25 


2 


5 


5 


4.75 


1 


3 


4 


5 


2 


1 


4 


1 


1.5 


1 


1 


5 


1 


1 


1 


1.25 



Table B.2: Wavelet: Real Daub-3, SNR: 10 dB 



S 


Real 


Imag. 


Angle 


Abs. 


1 


5 


5 


2.75 


1 


2 


5 


5 


5 


1 


3 


5 


5 


2 


1.25 


4 


1.5 


1 


1 


2.5 


5 


1 


1 


1 


1.5 



5 — sharply obvious with sharp boundaries. 

4 — obvious. 

3 — fair. 

2 — distinguishable from hackgroimd noise. 

1 — not distinguishable from background noise. 
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Table B.3: Wavelet: Complex Morlet, SNR: 10 dB 



s 


Real 


Imag. 


Angle 


Abs. 


1 


5 


5 


4.5 


1 


2 


5 


5 


4.5 


1 


3 


5 


5 


4.5 


1 


4 


4.25 


4.75 


1 


3.5 


5 


3.75 


3.75 


1 


2.75 



Table B.4: Wavelet: Real Morlet, SNR:10 dB 



S 


Real 


Imag. 


Angle 


Abs. 


1 


5 


5 


4.5 


4.5 


2 


5 


5 


4.5 


4.5 


3 


5 


5 


4.5 


4.75 


4 


4.75 


4.75 


1 


3.5 


5 


4.75 


3.75 


1 


3.75 



5 — sharply obvious with sharp boundaries. 

4 — obvious. 

3 — fciir. 

2 — distinguishable from background noise. 

1 — not distinguishable from backgroimd noise. 
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Table B.5: Wavelet: Complex Morlet, SNR:3 dB 



s 


Real 


Itnag. 


Angle 


Abs. 


1 


3.5 


3.25 


3.5 


4 


2 


4 


3.5 


4 


3.5 


3 


4.5 


4.5 


3.5 


3.5 


4 


2.5 


1.75 


1 


3 


5 


2 


2.25 


1 


3 



Table B.6: Wavelet: Real Morlet, SNR:3 dB 



S 


Real 


Imag. 


Angle 


Abs. 


1 


3 


1.75 


1 


3.75 


2 


4.5 


4 


2.75 


3.75 


3 


4.5 


4.5 


3 


4 


4 


2.75 


2.75 


1 


2 


5 


2.75 


2.75 


1 


1.5 



5 — sharply obvious with sharp boundaries. 

4 — obvious. 

3 — fair. 

2 — distinguishable from background noise. 

1 — not distinguishable from background noise. 
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Table B.7: Wavelet: Complex Daub-3, SNR:3 dB 



s 


Real 


Imag. 


Angle 


Abs. 


1 


2.8 


3.6 


2.2 


1 


2 


2.4 


3.6 


4.25 


2.8 


3 


3.4 


4.25 


1.4 


1 


4 


1.8 


1 


1 


1 


5 


1 


1 


1 


1 



Table B.8: Wavelet: Real Daub-3, SNR:3 dB 



S 


Real 


Imag. 


Angle 


Abs. 


1 


3.8 


1.4 


2.4 


1.6 


2 


3.6 


3.2 


4.4 


3.4 


3 


2.6 


3.2 


2.2 


1.2 


4 


1 


1 


1 


1 


5 


1 


1 


1 


1 



5 — sharply obvious with sharp boundaries. 

4 — obvious. 

3 — fair. 

2 — distinguishable from background noise. 

1 — not distinguishable from background noise. 



134 




APPENDIX C. WAVELET SURFACES OF 

OTHER SIGNALS 



This appendix includes the wavelet surfaces obtained by processing noise only, 
ASK, PSK and FSK signals. Figures 40 and 41 show no identifiable diamond structure 
pattern for the noise only case. Each pair of successive figures (Figure 42 and 43, 
Figure 44 and 45, and Figure 46 and 47) show a diamond structure residing only in 
one scale . 
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Re(S3) Re(S2) Re(S1) Re(CF) 



time delay 
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Figure 40. Wavelet surfaces obtained from noise only, scale 1, 2 and 3. 



136 




Re(S5) Re{S4) 



time delay 



Daub-3Noise only 




time 



Figure 41. Wavelet surfaces obtained from noise only, scale 4 and 5. 
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Figure 42. Wavelet surfaces obtained from ASK signal, scale 1, 2, and 3. 
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Figure 43. Wavelet surfaces obtained from ASK signal, scale 4 and 5. 
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Figure 44. Wavelet surfaces obtained from FSK signal, scale 1, 2 and 3. 
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Figure 45. Wavelet surfaces obtained from FSK signal, scale 4 and 5. 
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Figure 46. Wavelet surfaces obtained from PSK signal, scale 1, 2 and 3. 
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Figure 47. Wavelet surfaces obtained from PSK signal, scale 4 and 5. 
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